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This  dissertation  examines  the  electromagnetic  (EM)  properties  of  two  component 
composites.  These  composites  exhibit  varying  degrees  of  complexity,  ranging  from  those 
containing  single  cube  inclusions  to  those  consisting  of  hundreds  of  irregularly  shaped 
inclusions  in  a  disordered  arrangement.  Specifically,  time-domain  EM  modeling  is  used 
as  a  tool  to  calculate  effective  permittivity  as  well  as  EM  localization  behavior  for  these 
mixtures. 

For  single  inclusion  mixtures  with  high  permittivity  contrast  between  the  inclusion 
and  the  surrounding  matrix  material,  strong  deviation  from  established  mixing  theories  is 
found.  It  is  also  shown  that  the  orientation  of  the  inclusion  impacts  the  effective  electrical 
permittivity  of  the  composite.  Electric  fields  are  found  to  localize  on  the  edges  and  corners 
of  an  irregular  inclusion  independent  of  simulation  boundary  conditions. 

Increasing  the  complexity  of  the  studied  mixtures,  the  EM  properties  are  then  ana¬ 
lyzed  for  mixtures  containing  many  irregularly-  and  regularly-shaped  inclusion  crystals. 
A  strong  correlation  between  effective  permittivity  and  cross-sectional  area  is  found  for 
these  mixtures.  With  hundreds  of  inclusion  crystals  (with  either  irregular  or  cube  shapes), 
the  change  in  inclusion  shape  causes  negligible  differences  in  effective  permittivity. 

Electric  energy  density  localizes  on  edges  and  corners  of  inclusions  regardless  of  in¬ 
clusion  shape.  A  result  not  highlighted  by  single  inclusion  mixtures,  for  composites  con¬ 
taining  numerous  inclusions  significant  increases  in  energy  localization  are  observed  as 
the  EM  signal  travels  through  more  inclusions.  This  increase  in  hotspot  magnitude  is 
observed  for  complex  mixtures  with  hundreds  of  disordered  irregularly  shaped  inclusion 
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crystals.  The  same  behavior  is  also  observed  for  the  much  simpler  scenario  of  a  few  dozen 
cube  inclusions  distributed  within  discrete  parallel  planes.  As  a  result,  a  key  conclusion 
from  this  work  is  that  the  study  of  hotspots  does  not  require  the  use  of  geometrically 
complex  structures  with  hundreds  of  irregular  inclusions.  Rather,  much  simpler  arrange¬ 
ments  of  fewer  and  regularly  shaped  inclusions  can  be  analyzed  since  they  exhibit  similar 
energy  localization  behavior. 
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Chapter  1 


Introduction 

1.1  Overview 

Studies  of  the  analysis  and  detection  of  energetic  materials  [1]— [5]  are  of  substantial 
importance.  To  aid  in  this  research,  it  is  important  to  be  able  to  model  and  analyze 
component  mixtures.  To  understand  the  EM  characteristics  of  such  materials,  it  is  nec¬ 
essary  to  understand  the  coupling  of  EM  energy  into  random  mediums  with  specific 
structures  including  energetic  or  explosive  materials  and  organic  compounds.  Analysis  of 
the  concentration  of  EM  energy  is  significant  because  such  a  concentration  could  result 
in  increased  vapor  pressure  from  organic  composites,  making  it  easier  to  detect  these 
materials. 

Specifically,  it  takes  relatively  little  energy  to  release  energy  in  metastable  bonds  in 
a  material.  If  a  hotspot  is  introduced  into  an  energetic  material  it  can  heat  up.  The 
following  research  seeks  to  further  understand  the  behavior  of  such  localization,  as  well 
as  general  EM  properties  of  component  mixtures. 
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1.2  Motivations 


An  important  goal  of  this  research  is  to  further  explore  the  coupling  of  EM  energy  into 
disordered  and  random  mediums.  To  aid  in  this  understanding  the  objectives  of  this 
research  are  to  understand  if  these  hotspots  occur  at  certain  locations  and  what  condi¬ 
tions  lead  to  increased  hotspot  magnitudes.  Realistic  crystalline  composites  can  contain 
crystals  with  various  sizes,  and  these  crystals  can  be  both  regular  and  irregular  in  shape. 
So,  another  objective  of  this  research  is  to  determine  how  the  shape  and  distribution 
of  crystals  affects  how  energy  is  stored  within  a  mixture.  The  next  section  discusses  an 
important  type  of  mixture  that  is  especially  relevant  today,  namely  explosive  or  energetic 
materials. 

1.2.1  Explosive  Materials 

In  general  terms,  a  chemical  explosive  can  be  defined  as  “a  compound  or  a  mixture  of 
compounds  which,  when  subjected  to  heat,  impact,  friction,  or  shock,  undergoes  very 
rapid,  self-propagating,  heat-producing  decomposition”  [6].  EM  radiation  incident  on  an 
explosive  or  energetic  material  can  cause  a  release  of  energy  that  has  the  potential  to  be 
damaging  to  a  target  [6].  Energy  can  be  defined  as  either  potential,  which  measures  the 
ability  of  doing  work,  or  kinetic,  which  is  the  energy  of  motion  [6].  Explosives  such  as 
trinitrotoluene  (TNT)  have  chemical  bonds  that  are  relatively  unstable  [6].  So,  chemical 
potential  energy  exists  that  can  be  turned  into  kinetic  energy  under  the  right  conditions 
[6].  The  detonation  of  the  explosive  brings  stability  to  the  chemical  bonds,  leading  to 
kinetic  energy  achieved  through  blast  and  thermal  energy  release  [6].  Since  overall  there 
is  a  production  of  energy,  the  process  of  changing  the  bonds  of  a  chemical  system  to  a 
greater  amount  of  stability  is  an  exothermal  process  [6]. 
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High  explosives  (HEs)  will  detonate,  which  occurs  when  the  reaction  is  moving  through 
the  explosive  faster  than  the  speed  of  sound  in  the  explosive  that  has  not  reacted  [6]. 
Low-order  explosives  (LEs)  will  deflagrate  or  rapidly  burn  with  a  reaction  that  is  slower 
than  the  speed  of  sound  in  the  explosive  that  is  unreacted  [6].  Almost  all  commercially 
available  explosives,  not  including  black  powder,  are  HEs  [6]. 

Research  Department  Explosive  (RDX)  is  a  combination  of  a  crystalline  explosive  and 
a  binder,  which  is  often  a  plastic  material  [6].  High  Melt  Explosive  (HMX)  is  chemically 
similar  to  RDX  [6].  RDX  and  HMX  are  part  of  a  group  of  organic  nitrate  explosives 
called  nitramines  [7].  The  study  of  energetic  materials  is  an  important  movitation  of  this 
research.  HMX  is  of  widespread  importance,  being  the  “highest-energy  solid  explosive 
produced  on  a  large  scale  in  the  United  States”  [7].  In  addition,  the  Navy  uses  HMX  as 
a  plastic-bonded  explosive  ingredient  [7] . 

Improvised  Explosive  Devices  (IEDs)  are  of  significant  concern  in  places  such  as  Iraq 
[8] .  IEDs  in  Iraq  can  be  initiated  with  a  device  such  as  a  cell  phone  that  triggers  the  IED 
to  go  off  [8].  In  practice  there  are  several  methods  for  disarming  an  IED  including  using 
EM  pulses  [8].  As  a  method  for  either  disabling  or  destroying  an  IED,  a  high  energy  EM 
pulse  can  cause  the  trigger  of  an  IED  to  “fail  open”  and  not  induce  an  explosion  or  “fail 
closed”  and  detonate  the  IED  [8].  However,  very  cheap  countermeasures  can  be  used 
to  counter  attempts  at  destroying  or  disabling  IEDs  [8].  Clearly,  additional  methods 
are  required.  Since  EM  energy  can  be  used  to  destroy  or  disable  explosive  devices,  it 
is  important  to  gain  further  understanding  of  how  this  energy  behaves  within  random 
mixtures.  To  that  end,  an  important  motivation  of  this  research  is  the  need  to  study 
how  energy  is  stored  in  three-dimensional  random  mixtures,  including  those  mixtures 
containing  the  energetic  material  HMX.  HMX  is  a  type  of  polymer  bonded  explosive 
(PBX),  which  is  a  combination  of  a  polymeric  binder  with  explosive  crystals  [9].  The  EM 
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properties  of  HMX  crystals,  as  well  as  their  sizes  based  on  [9],  will  be  used  here  to  study 
realistic  crystalline  mixtures. 

A  motivation  in  studying  combinations  of  variously-sized  crystals  is  that  if  hotspots 
can  be  introduced  into  a  material  using  EM  signals,  the  crystal  phase  state  in  the  material 
can  be  altered.  This  has  various  implications  including  IED  destruction  or  disablement. 
The  energy  storage  capabilities  of  mixtures  will  also  be  investigated  with  studies  of  ef¬ 
fective  permittivity  of  composites. 

1.2.2  Broader  Impact 

The  research  described  in  this  dissertation  also  has  various  other  applications  that  extend 
beyond  energetic  materials.  One  such  application,  which  can  be  related  to  composites 
with  randomly  oriented  and  distributed  inclusions,  is  the  next  generation  of  aerospace- 
qualified  Carbon  Fiber  Reinforced  Composites  (CFRCs).  While  still  at  least  a  decade 
away,  these  composites  will  utilize  random  fibers  and  will  not  delaminate  (or  separate 
into  layers),  coming  with  a  greatly  reduced  cost  and  increased  reliability.  This  research 
can  also  be  related  to  several  other  applications  including  the  characterization  of  solid 
rocket  fuel,  the  demolition  of  either  visible  or  buried  structures,  remote  interrogation, 
and  demining. 

1.2.3  Technical  Approach 

In  this  research  CST  Microwave  Studio  [10]  will  provide  the  basis  for  studying  three- 
dimensional  mixtures.  CST  Microwave  Studio  is  a  software  tool  for  high  frequency  appli¬ 
cations  that  provides  3D  EM  simulations.  It  will  also  serve  as  a  useful  tool  for  studying 
the  permittivity  of  composites.  Specifically,  the  time-domain  solver  in  CST  Microwave 
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Studio,  utilizing  the  finite  integration  technique  (FIT),  will  be  utilized.  The  time-domain 
solver  allows  energy  distributions  to  be  studied  as  a  function  of  time. 

Utilizing  CST  Microwave  Studio,  materials  will  be  modeled  with  various  permittivi¬ 
ties,  shapes,  sizes,  and  orientations  to  show  how  these  characteristics  impact  the  permit¬ 
tivities  and  energy  storage  behavior  of  composite  mixtures.  The  crystals  that  are  part  of 
these  mixtures  will  be  given  material  properties  of  realistic  materials  including  HMX. 

1.3  Original  Contributions 

The  original  contributions  of  this  dissertation  are  summarized  here. 

1.  Two-component  composites  are  studied  involving  both  regularly-  and  irregularly- 
shaped  inclusions.  The  orientation  of  inclusions  is  shown  to  affect  the  effective 
electrical  permittivity  of  mixtures.  This  is  presented  in  Chapter  3. 

2.  Structures  are  analyzed  involving  many  scaled  and  rotated  inclusion  crystals  with 
both  regular  and  irregular  shapes.  These  complex  structures  are  simulated  in  CST 
Microwave  Studio  with  EM  analysis  of  energy  localization  behavior  and  effective 
permittivity.  This  is  presented  in  Chapter  4. 

3.  The  effective  permittivity  of  composite  mixtures  is  found  to  be  highly  correlated 
to  the  cross-sectional  area  of  the  inclusion  crystals  within  the  mixture.  This  is 
presented  in  Sections  3.6  and  4.6. 

4.  Effective  permittivity  is  compared  between  mixtures  involving  many  regularly-  and 
irregularly-shaped  inclusions.  With  hundreds  of  inclusions  the  effective  permittivity 
is  found  to  be  independent  of  inclusion  shape.  This  is  presented  in  Section  4.5. 
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5.  Both  disordered  and  more  organized  arrangements  of  inclusions  show  EM  localiza¬ 
tion  on  inclusion  edges  and  corners.  This  is  presented  in  Section  4.7. 

6.  Increased  hotspot  magnitude  is  found  in  composites  containing  many  inclusions  as 
the  EM  signal  travels  through  more  inclusions  (a  result  not  highlighted  by  single 
inclusion  mixtures).  This  is  presented  in  Section  4.7. 

1.4  Outline  of  Dissertation 

Chapter  2  of  this  dissertation  provides  a  literature  review  of  topics  related  to  effective 
permittivity  and  localization  of  EM  energy. 

Chapter  3  studies  two-component  composites  with  regularly-  and  irregularly-shaped 
inclusion  crystals.  This  chapter  focuses  on  the  relatively  simpler  scenarios  of  a  single 
irregular  inclusion  and  cube-shaped  inclusions  placed  on  an  ordered  grid.  With  high  per¬ 
mittivity  contrast  between  inclusion  and  matrix  material,  it  is  shown  that  the  orientation 
of  the  inclusion  affects  the  effective  relative  permittivity  of  the  mixture.  Fields  are  found 
to  localize  at  edges  and  corners  of  the  irregular  inclusion  independent  of  the  simulation 
boundary  conditions. 

Chapter  4  investigates  the  EM  properties  of  complex  mixtures  involving  many  crystals 
of  regular  and  irregular  shapes.  The  Bullet  physics  library  [11],  a  physics  engine  commonly 
used  for  computer  games,  is  used  to  create  randomized  combinations  of  many  crystals. 
This  automated  procedure  leads  to  a  higher  packing  of  crystals  than  can  be  obtained 
reasonably  with  manual  creation.  With  high  permittivity  contrast  between  inclusions 
and  matrix  material,  effective  permittivity  is  found  to  deviate  from  established  mixing 
theories.  With  hundreds  of  irregular  or  cube  inclusion  crystals  within  a  mixture,  the 
impact  of  individual  inclusion  shape  averages  out.  Surrounding  inclusions  are  found  to 
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impact  the  electric  energy  localization  within  the  mixtures,  with  higher  localization  seen 
after  the  EM  signal  has  traveled  through  more  inclusion  crystals. 

Chapter  5  provides  a  summary  of  the  research  presented  and  suggests  areas  for  future 
study. 

1.5  Published  Works 

The  work  outlined  in  this  dissertation  led  to  the  following  journal  papers. 

1.5.1  Journal  Papers 

1.  ©  2013  IEEE.  Reprinted,  with  permission,  from  A.  J.  Pickles  and  M.  B.  Steer, 
“Effective  permittivity  of  3-D  periodic  composites  with  regular  and  irregular  inclu¬ 
sions,”  IEEE  Access,  vol.  1,  pp.  523-536,  2013. 

2.  ©  2013  IEEE.  Reprinted,  with  permission,  from  A.  J.  Pickles,  I.  M.  Kilgore,  and  M. 
B.  Steer,  “Automated  creation  of  complex  three-dimensional  composite  mixtures 
for  use  in  electromagnetic  simulation,”  IEEE  Access,  vol.  1,  pp.  248-251,  2013. 

3.  ©  2013  IEEE.  Reprinted,  with  permission,  from  A.  J.  Pickles,  I.  M.  Kilgore,  and  M. 
B.  Steer,  “Electromagnetic  properties  of  disordered  three-dimensional  mixtures,” 
IEEE  Access,  submitted  for  publication,  September  2013. 
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Chapter  2 


Literature  Review 

2.1  Introduction 

The  main  topics  explored  in  the  following  chapters  are  effective  permittivity  and  EM 
energy  localization.  This  chapter  begins  by  introducing  the  concepts  of  permittivity  and 
effective  permittivity.  An  overview  is  provided  of  the  study  of  effective  permittivity  over 
the  past  few  centuries,  culminating  with  the  sophisticated  computer  software  used  today 
for  studying  effective  permittivity.  Finally,  this  chapter  presents  an  overview  of  research 
related  to  studies  of  effective  permittivity  and  energy  localization. 

2.2  Permittivity 

2.2.1  Introduction 

First,  this  section  defines  the  concept  of  permittivity  as  referenced  in  this  dissertation, 
beginning  with  an  overview  of  the  related  topics  of  electric  fields  and  dielectrics. 


2.2.2  Electric  Field,  Dielectrics,  and  Permittivity 

If  a  point  charge,  Q,  exists  and  experiences  a  force,  F,  due  to  additional  charges,  the 
electric  held,  E,  is  defined  as  [12] 

e=Fq-  P'1* 

From  Equation  2.1,  the  electric  held  direction  and  the  direction  of  the  force  exerted  by 
the  held  on  a  positively  charged  particle  are  the  same  [12].  The  attraction  of  opposite 
charges  and  the  repulsion  of  like  charges  means  the  electric  held  points  from  positive 
charges  and  to  negative  charges  [12]. 

When  an  electric  held  is  applied  to  a  group  of  materials  known  as  dielectrics  [13], 
which  are  insulators,  the  dielectric  can  be  polarized.  Dielectric  polarization  occurs  when 
charges  shift  from  their  original  positions  due  to  the  applied  electric  held  [13].  There  is  not 
a  how  of  charges  like  in  a  conductor,  but  rather  the  dielectric  polarization  causes  positive 
charges  to  move  towards  the  held  and  negative  charges  to  move  away  from  the  held 
[13].  This  newly  formed  separation  of  charges  in  the  dielectric  creates  an  internal  electric 
held  that  points  in  the  opposite  direction  to  the  applied  electric  held,  thus  reducing 
the  overall  held  in  the  dielectric  [13].  A  dielectric  material  has  energy  storing  capacity 
through  polarization,  such  as  the  material  between  a  capacitor’s  metallic  plates  [13] .  The 
surface  charge  of  the  capacitor  is  increased  when  an  electric  held  is  applied  because  the 
dielectric  becomes  polarized  [13]. 

The  term  “insulator”  generally  is  implied  to  have  a  low  amount  of  electrical  conduc¬ 
tion,  while  “dielectrics”  have  a  high  amount  of  polarizability,  measured  by  the  dielectric 
constant  [13].  The  dielectric  constant,  or  static  relative  permittivity,  is  the  relative  per¬ 
mittivity  exhibited  by  a  material  when  at  zero  frequency  [14].  However,  dielectric  constant 
could  refer  to  the  static  or  frequency-dependent  relative  permittivity  based  on  context 
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[14]- 

Permittivity,  e,  is  a  measure  of  a  material’s  ability  to  “permit”  an  electric  field  [15]. 
The  permittivity  in  vacuum  is  eo  [14].  The  relative  permittivity,  er,  as  a  function  of 
angular  frequency,  u,  is  defined  as  [14] 


From  Equation  2.2,  the  relative  permittivity  of  a  material  is  a  ratio  of  how  much  electric 
energy  is  stored  in  the  material  when  a  voltage  is  applied,  e(oj),  relative  to  how  much 
is  stored  in  a  vacuum,  eo  [14].  Relative  permittivity  is  generally  complex,  and  the  imag¬ 
inary  part  of  the  permittivity  causes  EM  attenuation  inside  the  medium  [14].  Relative 
permittivity  for  materials  with  loss  is  defined  as  [14] 


In  Equation  2.3,  ereai  is  the  real  part  of  the  relative  permittivity  and  a  is  the  conductivity 
[14].  The  loss  tangent  is  defined  as  the  ratio  of  the  imaginary  over  the  real  part  of  the 
permittivity  [16].  Having  defined  the  concept  of  permittivity,  the  next  section  provides 
an  overview  of  how  the  study  of  permittivity  has  evolved. 


2.3  Effective  Permittivity 

2.3.1  Introduction 

While  permittivity  is  a  useful  parameter  to  describe  the  electric  energy  storage  ability 
of  individual  materials,  as  summarized  in  the  previous  section,  an  additional  concept  is 
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necessary  to  describe  the  overall  permittivity  of  a  mixture  that  contains  more  than  one 
material.  This  concept  is  described  by  effective  permittivity.  Exploration  of  the  effective 
permittivity  of  a  composite  medium  is  a  major  focus  of  the  research  presented  in  the 
following  chapters  of  this  dissertation. 

This  section  provides  an  overview  of  how  the  study  of  effective  permittivity  has 
evolved  over  the  past  few  centuries,  as  discussed  in  detail  in  [17]  and  summarized  here. 
These  stages  include  [17]  the  work  of  James  Clerk  Maxwell,  the  Maxwell  Garnett  mixing 
law,  Bruggeman’s  work  on  the  effective  medium,  bounding  methods,  and  finally  compu¬ 
tation  of  electromagnetics  with  the  use  of  computers.  The  use  of  sophisticated  computer 
software  is  especially  relevant  here  since  it  is  used  for  the  research  presented  in  the  fol¬ 
lowing  chapters. 

2.3.2  Early  Work 

In  the  late  19th  century,  James  Clerk  Maxwell  was  among  the  first  to  develop  formulas 
for  the  effective  permittivity  of  heterostructures  [17].  Maxwell  [18]  studied  the  variation 
in  potential  for  the  simple  scenario  of  a  spherical  inclusion  of  one  material  inside  another 
material  [17].  Soon  after,  formulae  were  developed  to  describe  approximately  the  proper¬ 
ties  of  the  effective  medium  with  knowledge  of  the  materials  that  made  up  the  medium 
[17].  Gladstone  and  Dale  [19]  utilized  an  approximation  formula  where  the  permittivity 
of  the  mixture  was  linearly  proportional  to  the  volume  taken  up  by  and  permittivity  of 
the  components  [17]. 

Analytic  approaches,  called  mixing  laws,  to  determining  effective  permittivity  rely 
on  a  volume  average  of  individual  permittivities  [17].  While  the  volume  fractions  are 
important  parameters  of  mixing  laws,  internal  sizes  of  the  inclusions  are  not  considered 
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[17].  The  accuracy  and  applicability  of  these  laws  has  been  questioned  as  they  have 
been  limited  to  relatively  simple  scenarios  that  are  dilute  and  geometrically  simple  [17], 
[20].  This  dissertation  furthur  studies  the  applicability  of  mixing  laws  by  comparing  the 
simulated  results  to  those  predicted  by  the  famous  Maxwell  Garnett  mixing  law  [21]. 

2.3.3  Effective  Medium  Theories 

Additional  understanding  of  the  properties  of  mixtures  came  with  the  development  of 
the  effective  medium  approximation  (EMA),  significantly  discussed  by  Bruggeman  [22] 
in  1935  [17].  Overall,  “EMA  approaches  assume  that  each  constituent  is  surrounded  by 
the  same  effective  medium.  It  assumes  that  the  local  electric  and  magnetic  fields  are  the 
same  in  the  volume  occupied  by  each  component  in  the  composite.  The  basic  understand¬ 
ing  appears  to  be  that  the  fictitious  replacement  (effective)  medium  is  macroscopically 
equivalent  to  the  real  medium  in  its  dielectric  effects.  The  analysis  is  done  in  the  approx¬ 
imation  of  noninteracting  inclusions  (each  inclusion  is  subject  to  the  same  mean  held, 
unperturbed  by  the  presense  of  other  inclusions)”  [17].  EMA  methods  are  applicable  at 
low  volume  fractions  of  inclusions  with  the  assumptions  that  groupings  or  overlapping  of 
inclusions  do  not  occur  [17].  Since  in  the  effective  medium  the  energy  density  is  seen  to  be 
uniform  throughout  the  entire  volume,  EMA  methods  break  down  at  higher  volume  frac¬ 
tions  of  inclusions  since  they  do  not  take  into  account  the  interactions  among  the  various 
components  of  a  mixture  [17].  So,  EMA  approaches  do  not  provide  an  all-encompassing 
complete  solution  for  the  effective  permittivity  of  two-component  mixtures  [17]. 

In  regards  to  frequency,  effective  medium  theories  (EMTs)  are  valid  in  the  quasistatic 
limit,  meaning  that  the  wavelength  is  much  larger  than  the  linear  dimensions  of  the 
individual  inclusions  that  make  up  the  mixture  [17].  So,  although  the  wave  cannot  ‘see’ 
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the  inclusions,  the  effective  permittivity  is  still  impacted  by  the  inclusions  [17]. 

2.3.4  Bounds  on  Permittivity 

Another  important  advance  in  the  theory  of  effective  permittivity  came  with  the  estab¬ 
lishment  of  rigorous  bounds  on  the  effective  permittivity  of  two-component  mixtures  [17]. 
Absolute  bounds  for  effective  permittivity  were  developed  by  Wiener  [23],  occurring  for 
situations  of  inclusion  layers  perpendicular  or  parallel  to  the  applied  held  [17]. 

2.3.5  Electromagnetics  Aided  by  Computers 

The  development  of  computers  has  advanced  EM  study.  Known  as  computational  elec¬ 
tromagnetics  (CE),  this  area  connects  electromagnetism  with  materials  science  and  com¬ 
puter  science  [17].  Results  taken  from  computer  simulations  can  be  used  to  confirm  or 
disprove  EMTs,  and  suggest  new  areas  that  deserve  further  research  [17].  Numerous  CE 
methods  [17]  have  been  used  in  effective  permittivity  investigations.  These  include  Fourier 
expansion  [24],  [25],  Finite  Difference  Time  Domain  (FDTD)  [26],  [27],  [28],  finite  element 
(FE)  [29],  [30],  [31],  the  transfer  matrix  method  [32],  and  boundary-integral  equations 
(BIE)  [33],  [34],  [35].  In  a  basic  sense,  these  methods  involve  a  combination  of  Maxwell’s 
equations  with  a  set  of  boundary  conditions  [17].  Sophisticated  computer  software,  such 
as  CST  Microwave  Studio,  can  be  used  to  gain  greater  understanding  of  the  EM  behavior 
of  complex  component  mixtures. 

2.3.6  Summary 

As  outlined  in  this  section,  the  development  of  mixing  laws  and  EMAs  were  important 
milestones  in  the  study  of  effective  permittivity.  The  results  of  the  simulated  effective 
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permittivities  in  the  following  chapters  will  be  compared  to  the  Maxwell  Garnett  and 
Bruggeman  theories  to  further  explore  their  applicability. 

2.4  Research  on  Composite  Mixtures 

2.4.1  Introduction 

The  previous  section  covered  previous  research  related  to  mixing  laws  and  EMAs.  This 
section  explores  past  and  current  research  as  it  relates  to  the  subject  of  effective  permit¬ 
tivity  and  EM  energy  localization  in  composite  mixtures. 

2.4.2  Effective  Permittivity  Research 

Studies  of  effective  permittivity  have  included  mixtures  organized  in  both  two-  and  three- 
dimensions.  Research  has  been  performed  looking  at  component  mixtures  where  the 
inclusions  have  regular  and  symmetric  shapes.  Karkkainen,  Sihvola,  and  Nikoskinen  [36] 
studied  two-dimensional  mixtures  with  inclusions  in  the  shape  of  discs.  The  authors 
utilized  mixtures  contained  within  a  TEM  (Transverse  ElectroMagnetic)  waveguide  and 
the  inclusions  were  allowed  to  overlap  with  each  other  [36].  They  considered  only  one 
permittivity  contrast  between  matrix  and  inclusion  material  and  did  not  find  a  model 
that  predicted  the  simulated  results  over  all  volume  fractions  [36]. 

For  a  TEM  wave,  the  electric  held,  magnetic  held,  and  direction  of  propagation  are 
all  orthogonal  to  each  other.  An  example  of  the  axis  definitions  for  a  TEM  environment 
is  given  in  Figure  2.1.  In  Figure  2.1,  the  electric  held  is  polarized  in  the  y  direction,  the 
magnetic  held  is  polarized  in  the  z  direction,  and  propagation  is  in  the  x  direction. 

In  another  paper,  Karkkainen,  Sihvola,  and  Nikoskinen  [37]  studied  three-dimensional 
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Figure  2.1:  Axis  definitions  for  a  TEM  wave. 


mixtures  of  spheres  and  found  different  results  for  effective  permittivity  based  on  whether 
or  not  the  spheres  were  allowed  to  form  clusters.  A  mixing  model  was  developed  to  fit 
the  simulated  results  based  on  least  squares  fitting  [37] . 

Wu,  Chen,  and  Liu  [38]  observed  differences  in  effective  permittivity  with  spherical 
inclusions  compared  to  inclusions  shaped  like  a  cube.  Wu  et  al.  [39]  also  studied  spheri¬ 
cal  inclusions  in  three  dimensions.  The  authors  found  good  agreement  with  conventional 
mixing  theories  when  there  was  low  contrast  between  the  sphere  and  surrounding  mate¬ 
rial. 

Luo,  Liu,  and  Pan  [40]  analyzed  three-dimensional  mixtures  with  spherical  inclusions 
as  well  as  spherical  inclusions  with  an  outer  membrane.  The  authors  studied  how  features 
such  as  the  membrane  thickness  had  a  strong  impact  on  effective  permittivity  and  found 
differences  with  different  inclusion  shapes  including  cubes  and  cylinders  [40].  Wu  et  al. 
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[41]  studied  differences  in  effective  permittivity  with  different  spacial  arrangements  of 
spherical  inclusions  as  well  as  a  cylindrical  inclusion. 

Cheng  et  al.  [42]  analyzed  three-dimensional  mixtures  where  cubes  of  two  different 
permittivities  were  organized  in  a  random  arrangement  to  create  a  larger  cube-shaped 
mixture.  Since  the  larger  created  structure  was  anisotropic  in  three  dimensions,  its  effec¬ 
tive  permittivity  was  dependent  on  orientation  [42],  Chen  et  al.  [43]  created  a  structure 
of  cube  cells,  each  cell  with  a  sphere  inside.  The  sphere  was  then  given  a  randomly  vary¬ 
ing  radius  to  represent  different  sized  inclusions  within  the  entire  mixture  [43].  Finally, 
Jylha  and  Sihvola  [44]  studied  three-dimensional  mixtures  of  equal-sized  spheres.  The 
spheres  were  allowed  to  exist  in  random  positions  within  the  mixture  [44],  The  authors 
used  a  scalar  fitting  parameter  to  match  their  simulated  results  and  compared  the  sim¬ 
ulated  effective  permittivity  to  conventional  mixing  theories  [44],  In  [42],  [43],  and  [44], 
the  authors  used  an  averaging  of  the  effective  permittivity  determined  in  each  of  three 
orthogonal  directions  to  obtain  a  single  effective  permittivity  value. 

2.4.3  Energy  Localization  Research 

In  addition  to  effective  permittivity,  prior  research  has  also  investigated  EM  field  localiza¬ 
tion  behavior.  Mejdonbi  and  Brosseau  analyzed  regular  polygon  [45]  and  fractal  shaped 
[45],  [46]  inclusions  in  two  dimensions  using  FDTD  modeling.  Specifically  studying  an 
inclusion  with  a  fractal  shape,  the  authors  found  electric  field  hotspots  at  the  inclusion 
corners  [45].  The  authors  noted  that  these  hotspots  showed  a  lack  of  uniformity  in  field 
distributions  that  would  impact  effective  permittivity  [45].  In  addition,  EMA  methods 
based  on  uniform  field  distrubitions  (as  discussed  in  Section  2.3.3)  do  not  consider  these 
hotspots  and  so  are  limited  in  their  ability  to  describe  the  EM  properties  of  composites 
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[46], 

Using  the  FE  method,  Mejdoubi  and  Brosseau  [47]  analyzed  electric  field  localization 
on  fractal  shaped  inclusions.  These  hotspots  were  found  on  the  inclusion  corners  and 
edges,  showing  that  the  fields  were  not  uniform  throughout  the  composite  [47].  Also 
using  the  FE  method,  Cheng  et  al.  [42]  studied  three-dimensional  arragements  of  cube 
inclusions.  Each  cube  inclusion  within  a  10  x  10  x  10  arrangement  was  randomly  allowed 
to  have  one  of  two  permittivity  values,  creating  a  nonuniform  structure  [42] .  The  authors 
found  electric  field  hotspots  at  various  locations  within  the  stucture  [42], 

Electric  field  localization  can  be  studied  in  either  the  time-  or  frequency-domains. 
In  this  dissertation  the  time-domain  FIT  is  specifically  chosen  because  it  allows  field 
localization  behavior  to  be  studied  as  a  function  of  time.  Probes  can  be  placed  within 
the  structure  at  the  specific  locations  of  hotspots  and  the  electric  fields  at  those  points 
can  be  measured  throughout  the  simulation.  The  result  is  a  record  of  varying  field  be¬ 
havior  in  time  that  cannot  be  seen  with  steady-state  frequency  domain  analyses.  In  this 
dissertation,  field  probes  will  allow  localization  behavior  to  be  studied  for  increasingly 
complicated  arrangements  of  inclusions  to  show  how  energy  localizes  in  complex  mixtures. 

2.4.4  Summary 

Previous  researchers  have  analyzed  effective  permittivity  and  energy  localization  in  com¬ 
posite  mixtures.  This  research  used  symmetric  and  regularly  shaped  inclusions,  and  ex¬ 
panding  this  research  to  inclusions  of  an  irregular  shape  is  an  important  component  of  the 
following  chapters.  Realistic  composites  such  as  those  containing  energetic  materials  are 
often  not  made  up  of  inclusions  that  are  symmetric  and  uniform  in  shape.  The  differences 
in  effective  permittivity  with  irregularly  shaped  inclusions  presents  a  scenario  with  more 
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realistic  structures  than  those  created  only  with  spheres.  In  addition,  the  limitations  of 
the  classical  mixing  theories  can  be  further  understood  by  increasing  the  irregularity  in 
inclusion  shape. 

2.5  Summary 

This  chapter  reviewed  research  efforts  related  to  effective  permittivity  and  energy  local¬ 
ization,  which  are  continued  with  the  following  chapters.  Effective  permittivity  and  energy 
localization  describe  the  behavior  of  EM  energy  within  a  medium.  With  an  increase  in 
the  complexity  of  the  component  mixtures  analyzed,  the  study  of  effective  permittivity 
and  EM  localization  behavior  will  bring  a  greater  understanding  of  EM  energy  propagat¬ 
ing  within  more  complex  component  mixtures.  The  next  chapter,  Chapter  3,  describes 
studies  of  mixtures  with  single  irregular  inclusions,  and  for  comparison  purposes,  ordered 
cube  inclusions  on  a  grid.  The  complexity  of  the  analyzed  mixtures  is  increased  in  Chap¬ 
ter  4  with  the  introduction  of  composites  containing  many  disordered  irregular  and  cube 
inclusions. 
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Chapter  3 


Electromagnetic  Properties  of 
Mixtures  with  Single  and  Regularly 
Ordered  Inclusions 

3.1  Introduction 

This  chapter  explores  the  effective  temporal  surface-illuminated  properties  of  two-comp¬ 
onent  composites  consisting  of  inclusions  of  regularly-  and  irregularly-shaped  crystals 
in  a  matrix.  Time-domain  EM  modeling  using  the  FIT  is  used  to  calculate  scattering 
(S-)  parameters,  and  from  these  the  effective  relative  permittivities  are  calculated.  It 
is  shown  that  the  orientation  of  inclusions  with  high  permittivity  contrast  affects  the 
effective  electrical  permittivity  of  a  composite  mixture.  For  both  low  and  high  contrast 
inclusions,  fields  localize  at  edges  and  corners  of  the  irregular  inclusion  in  a  manner  not 
dependent  on  boundary  conditions  used  in  simulation. 

The  detection  and  neutralization  of  energetic  materials  using  stand-off  techniques  is 
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of  significant  and  growing  importance  [1]— [5] .  Most  solid  energetic  materials  are  random 
mixtures  comprising  irregularly-shaped  explosive  crystals  embedded  in  a  passive  com¬ 
posite  or  binder.  In  part,  the  energetic  material  can  be  characterized  by  its  steady-state 
bulk  electrical  permittivity.  However,  with  standoff  EM  probing,  the  probe  will  interact 
with  the  surface  and  near-surface  structure  of  the  energetic  material  rather  than  the  ma¬ 
terial  bulk.  Also  the  probe  will  be  transient  with  the  EM  waveform  appearing  as  a  radio 
frequency  (RF)  pulse  as  the  probe  scans  the  energetic  material.  This  chapter  explores 
the  electrical  properties  and  energy  localization  of  irregular  materials  excited  by  a  tem¬ 
poral  EM  pulse.  The  electrical  properties  determined  are  compared  with  the  standard 
mixing  theories  based  on  extreme  mixing  situations  and  using  steady-state  considera¬ 
tions  alone.  Time-domain  monitors  are  also  used  to  analyze  the  near-maximum  temporal 
energy  distributions  for  these  irregular  crystals. 

Synthetic  composites  comprising  irregularly-shaped  crystals  in  a  matrix  are  examined 
here  with  a  view  to  characterizing  their  effective  EM  properties  and  so  establish  a  base¬ 
line  for  remote  detection  and  material  characterization.  In  the  past,  studies  of  effective 
medium  properties  have  used  time-  and  frequency-domain  EM  simulation  methods  with 
various  abstractions  and  simplifications  made  to  render  the  simulation  computationally 
tractable.  The  most  common  simplification  is  to  consider  regularly-shaped  inclusions  and 
sometimes  a  2-D  projected  structure  is  analyzed.  In  the  time  domain,  the  effective  per¬ 
mittivity  has  been  calculated  in  2-D  for  random  mixtures  using  FDTD  simulation  [36] 
(where  inclusions  are  disks),  [45]  (disks,  regular  polygons  and  fractals),  and  3-D  random 
mixtures  using  FDTD  [38]  (spheres,  crosses  and  regular  polyhedrons),  and  [39]  (spheres 
and  a  complex  yet  symmetric  structure).  With  time-domain  EM  analysis,  regular  shapes, 
such  as  cubes,  spheres,  and  regular  polygons  have  been  studied  because  these  enabled 
simplified  gridding  using  rectangular  meshing  (and  analytic  projections  onto  a  rectangu- 
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lar  mesh  in  the  case  of  spheres). 

Frequency-domain  EM  analysis  of  2-D  random  mixtures  of  regularly- shaped  inclusions 
has  been  studied  using  the  finite  element  method  [47]  (polygons  and  fractals);  3-D  random 
mixtures  of  regularly-shaped  inclusions  using  a  frequency-domain  finite  difference  method 
[37]  (spheres),  [40]  (spheres,  crosses,  cylinders,  and  polyhedrons),  and  [41]  (spheres  and 
cylinders);  and  3-D  random  mixtures  using  a  frequency-domain  finite  element  method 
[42]  (polyhedrons),  and  [43],  [44]  (spheres).  Such  steady-state  analyses  are  limited  as  they, 
for  example,  do  not  capture  temporal  localization  of  energy. 

In  remote  probing,  a  region  is  scanned  and  a  steady-state  is  not  necessarily  obtained 
as  the  applied  signal  is  effectively  a  radio  frequency  (RF)  pulse.  The  significant  internal 
reflections  of  a  crystal-based  composite  illuminated  by  an  RF  pulse  result  in  what  is 
referred  to  as  a  long-tail  response.  This  characteristic  is  similar  to  that  observed  with  RF 
bandpass  filters  where  multiple  internal  reflections  in  response  to  an  RF  pulse  cannot  be 
predicted  from  steady-state  observations  [48]. 

A  composite  of  two  materials  will  typically  be  anisotropic  due  to  the  formation,  fab¬ 
rication,  or  manufacturing  process  favoring  particular  crystal  orientations.  When  the 
crystals  are  much  smaller  than  a  wavelength,  in  which  case  anisotropy  results  from  the 
combination  of  shapes,  the  standard  characterization  procedure  is  to  determine  an  ef¬ 
fective  permittivity  of  the  mixture  by  averaging  the  extracted  effective  permittivities  of 
each  of  three  orthogonal  directions  [42],  [43],  [44], 

This  chapter  studies  two  materials  comprising  crystals  in  a  lossless  matrix.  The  first 
is  an  explosive  material  consisting  of  HMX  crystals  [9]  and  the  other  is  a  material  with 
crystals  having  a  higher  relative  permittivity.  For  the  high  permittivity  material,  both 
regular  and  irregular  shapes  are  studied.  For  both  materials,  realistic  crystal  shapes  (as 
opposed  to  spheres,  cubes,  etc.)  are  considered.  As  such,  the  inclusions  do  not  necessar- 
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ily  have  symmetric  shapes.  Effective  properties  are  extracted  for  the  composites  using 
forward  and  reverse  propagation,  and  for  various  crystal  orientations.  The  overall  goals 
of  this  research  are  to  determine  under  what  situations  irregularly  shaped  inclusions  be¬ 
have  similarly  to  spheres  (through  comparison  with  the  classical  mixing  rules  based  on 
spheres)  and  cubes  in  calculations  of  effective  permittivity,  and  to  analyze  how  rotation 
of  the  inclusion  influences  the  effective  permittivity  of  composites.  EM  modeling  of  mix¬ 
tures  with  irregularly  shaped  inclusions  is  computationally  intensive  and  knowing  when 
simpler  geometries  can  be  used  is  advantageous. 

3.2  Three-Dimensional  Mixing  Rules 

Several  conventional  mixing  theories  exist  that  predict  the  effective  properties  of  two- 
component  mixtures.  One  of  these  is  the  Maxwell  Garnett  mixing  theory,  which  provides 
the  effective  permittivity  of  a  composite  as  [21] 


Off  —  e2  +  3ge2 — — - t - r-  (3-1) 

G  +  2e2  —  q(e  i  —  62) 

Here  is  the  permittivity  of  spherical  inclusions,  62  is  the  permittivity  of  the  surrounding 
material,  and  q  is  the  volume  fraction,  i.e.  filling  factor,  of  the  inclusions.  The  inclusion 
size  here  is  much  smaller  than  the  EM  wavelength.  A  few  assumptions  of  the  Maxwell 
Garnett  mixing  theory  are  that  the  inclusions  are  spherical  in  shape,  the  spheres  do  not 
touch  each  other,  and  the  radius  of  an  individual  spherical  inclusion  is  much  smaller  than 
the  distances  between  the  spheres  [42],  That  is,  the  mixing  rule  can  only  be  properly  used 
with  composites  that  have  a  low  proportion  of  inclusions  whereas  the  crystals  in  energetic 
materials  can  have  a  volume  fraction  of  up  to  95%. 
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The  second  classical  mixing  theory  is  known  as  the  Bruggeman  rule  [22] 

(3. 

As  with  the  Maxwell  Garnett  mixing  theory,  the  Bruggeman  mixing  rule  is  derived  with 
the  assumption  of  widely  separated  spherical  inclusions.  More  complicated  structures, 
such  as  irregularly-shaped  crystals,  are  expected  to  provide  different  results  [42], 

As  described  in  [57] ,  an  important  assumption  of  the  Bruggeman  rule  is  that  there  is  no 
separation  between  the  matrix  material  and  the  inclusions.  Rather,  the  mixture  is  treated 
as  a  homogenized  medium  over  which  polarizations  are  evaluated.  Thus  the  Bruggeman 
rule  treats  the  inclusions  and  surrounding  material  as  being  symmetric.  The  Bruggeman 
rule,  Equation  3.2,  balances  the  components  in  relation  to  the  effective  medium,  weighting 
each  component  by  q  for  inclusions  and  1— q  for  the  surroundings.  Conversely,  the  Maxwell 
Garnett  theory  does  not  utilize  such  symmetry. 

Mixing  rules  developed  for  composites  are  based  on  the  two  classical  mixing  rules 
(Maxwell  Garnett  and  Bruggeman)  but  none  cover  all  types  of  inclusions  [42],  In  this 
chapter  the  applicability  and  limitations  of  the  classic  Maxwell  Garnett  and  Bruggeman 
rules  are  further  developed. 
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3.3  Simulation  Method  and  Technique 

The  effective  permittivity  of  a  sample  of  finite  thickness  can  be  determined  using  mea¬ 
sured  or  simulated  scattering  (S-)  parameters  as  described  in  [49],  [50],  [51],  [52],  The 
time-domain  FIT  solver  with  hexahedral  meshing  implemented  in  the  commercial  EM 
package  GST  Microwave  Studio  is  utilized  for  these  simulations.  This  package  allows 
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complex  3D  shapes  to  be  modeled  and  simulated.  FIT  is  equivalent  computationally  to 
FDTD  except  that  the  time-domain  Maxwell’s  Equations  are  discretized  in  integral  form 
rather  than  differential  form  [52]  and  gridding  is  generalized. 

The  adaptive  meshing  algorithm  used  refines  the  mesh  until  convergence  of  electrical 
characterization  is  obtained.  Refinement  of  the  time-domain  mesh  by  an  additional  step 
led  to  an  effective  permittivity  change  of  less  than  0.5%  (indicating  that  acceptable 
accuracy  had  been  reached).  The  number  of  mesh  cells  for  each  time-domain  simulation 
run  was  at  least  250,000.  Computation  was  performed  using  an  80  core  cluster  with 
160  GB  of  RAM  and  clocking  at  2.66  GHz.  Where  required,  time  results  were  Fourier 
transformed  to  yield  a  frequency-domain  characterization  from  an  applied  wideband  EM 
pulse. 

With  the  RF  component  having  angular  frequency  u,  the  free  space  wavenumber 
k  =  lo/c  where  c  is  the  speed  of  light.  With  excitation  at  Port  1  (the  incidence  port),  the 
wavenumber,  the  thickness  d  of  the  structure  in  the  direction  of  propagation,  and  the 
calculated  A-parameters  were  used  to  derive  the  refractive  index  n  of  a  sample  surrounded 
by  vacuum  using  [51],  [52] 


n  —  ± 


1 

— arccos 
kd 


1 

25^ 


(1  -  S2n  +  si) 


(3.3) 


where  m  is  an  integer  indicating  multiple  possible  solutions.  In  Equation  3.3,  n  is  negative 
if  the  permittivity  and  permeability  are  both  negative.  The  relative  wave  impedance  z 
of  the  sample  is  [51],  [52] 


(1  +  Sn)2  -  SI 

(i  -  s„)’ -  s2y 


(3.4) 


The  S-parameters  in  Equation  3.3  and  Equation  3.4  are  referenced  to  the  impedance  of 
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free  space,  rj.  Information  about  the  material  can  be  used  to  remove  the  ambiguities  in 
Equation  3.3  and  Equation  3.4.  When  d  is  large  compared  to  the  wavelength  A  in  the 
medium,  the  branches  of  the  arccosine  function  determined  by  m  in  Equation  3.3  can  be 
close  together.  This  can  make  hireling  a  unique  solution  for  n  difficult.  However,  in  the 
following  simulations  d  is  less  than  A/4  so  that  there  is  a  unique  result  for  n.  As  a  result, 
the  default  branch  of  m  =  0  and  positive  n  is  chosen  in  Equation  3.3.  A  passive  material 
means  that  9fJ(z)  is  positive,  so  that  the  positive  branch  is  chosen  for  Equation  3.4.  Now, 
with  unique  solutions  for  n  and  z,  the  effective  permittivity  is  [51],  [52] 

eeff  =  (3.5) 

Here  n  can  be  well  defined  by  a  material  that  supports  one  propagating  mode  at 
a  certain  frequency  [51],  which  is  true  here.  For  the  frequency  range  studied  of  up  to 
50  GHz,  there  is  only  one  propagating  mode,  TEM,  all  other  modes  are  evanescent,  and 
Equations  3.3  and  3.4  use  the  transmission  and  reflection  coefficients  for  the  mode  that 
is  propagating.  Thus,  the  bulk  effective  permittivity  relies  on  an  abstraction  to  utilize 
the  single  propagating  mode  rather  than  on  evanescent  modes  in  its  derivation.  The 
results  that  follow  for  the  effective  permittivities  of  mixtures  use  the  procedure  above  to 
determine  eeff  unambiguously. 

3.4  Modeling  of  Irregular  Crystals 

Several  different  types  of  irregular  inclusions  are  analyzed  here  in  the  following  subsec¬ 
tions. 
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3.4.1  Modeling  of  Irregular  Crystals  in  CST 

An  important  goal  of  this  research  is  modeling  inclusions  of  a  size  and  shape  represen¬ 
tative  of  materials  comprising  irregularly-shaped  crystal  inclusions.  For  reference,  the 
crystal  size  of  the  irregular  energetic  material  HMX  is  typically  0.15  mm  [9].  This  was 
modeled  in  CST  as  the  irregularly-shaped  crystal  shown  in  Figure  3.1.  This  becomes  an 
inclusion  in  a  matrix  modeled  here  as  having  a  relative  permittivity  of  1.  In  Figure  3.1, 
the  irregular  crystal  spans  0.15  mm  in  the  z  direction  with  a  maximum  dimension  of 
0.178  mm  and  0.192  mm  in  the  x  and  y  directions,  respectively.  The  total  volume  of  the 
irregular  crystal  is  0.00351  nun3.  The  outer  box  shown  is  a  cube  with  a  side  length  of 
0.26  mm,  giving  an  inclusion  volume  fraction  of  q  =  0.20  for  Figure  3.1.  For  propagation 
in  free  space,  a  unit  cell  side  length  of  0.26  mm  corresponds  to  approximately  0.00087A 
(A  is  the  wavelength  measured  in  mm)  at  the  studied  frequency  of  1  GHz.  The  size  of 
the  irregular  crystal  is  kept  constant  and  the  outer  box  side  length  is  changed  to  obtain 
the  various  volume  fractions. 

EM  propagation  through  the  material  is  modeled  using  CST  Microwave  Studio  with 
the  crystal  placed  in  a  parallel  plate  transverse  EM  (TEM)  environment  [45].  With  the 
electric  held  polarized  in  the  x  direction,  the  TEM  environment  is  established  by  setting 
the  opposite  pair  of  normal  faces  in  the  x  direction  as  electric  walls,  the  opposite  pair  of 
normal  faces  in  the  y  direction  as  magnetic  walls,  and  the  final  opposite  pair  of  normal 
faces  in  the  z  direction  using  open  boundary  conditions.  In  simulation  the  open  faces 
are  modeled  as  perfectly  matched  layers  (PMLs)  that  eliminate  reflections  back  into  the 
TEM  environment.  Measurement  ports  are  set  up  at  the  open  faces  (Port  1  and  Port  2), 
defining  the  axis  of  EM  propagation  in  the  z  direction,  see  Figure  3.2. 

In  Figure  3.2,  the  ports  are  separated  from  the  mixture,  providing  an  averaging  effect 


26 


Figure  3.1:  An  irregularly-shaped  inclusion  inside  an  outer  box.  The  inclusion  volume 
fraction  here  is  20%.  The  crystal  is  within  a  cube  with  a  side  length  of  0.26  mm.  For 
comparison,  with  er  =  28,  the  wavelength  at  50  GHz  is  1.13  mm. 
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Figure  3.2:  The  irregular  crystal  of  Figure  3.1  shown  with  waveguide  ports  (Port  1  and 
Port  2)  set  along  the  propagation  axis,  the  z  axis.  In  turn  the  x  and  y  propagation 
directions  were  also  considered  with  the  crystal  structure  fixed  in  position  and  boundary 
conditions  changed  appropriately. 


similar  to  what  happens  with  stand-off  remote  sensing.  Specifically,  the  spacing  from  the 
waveguide  ports  to  the  structure  is  A/4  of  the  maximum  frequency  component  (50  GHz) 
of  the  applied  RF  pulse.  Propagation  in  free  space  results  in  a  separation  distance  of  1.5 
mm  for  both  ports  in  Figure  3.2.  Measurements  at  the  ports  shown  in  Figure  3.2  were 
referred  to  reference  planes  on  the  surface  of  the  cube  also  shown  in  Figure  3.2.  This 
yielded  the  S'-parameters  used  with  Equation  3.3  and  Equation  3.4  to  determine  the 
effective  permittivity  of  the  mixture  calculated  using  a  Gaussian  pulse  excitation  signal, 
Fourier  analyzed  temporal  response,  and  Equation  3.5. 

One  of  the  issues  in  EM  modeling  is  that  the  size  of  the  structure  that  can  be  simulated 
is  limited  and  it  is  necessary  to  establish  simulation  boundaries  close  to  the  structure 
being  investigated.  Thus  the  possible  impact  that  the  boundary  conditions  have  on  the 
extracted  results  is  a  concern.  The  TEM  waveguide  environment  is  necessary  to  extract 
the  effective  permittivity  and  the  environment  (i.e.,  the  TEM  set  of  boundary  conditions) 
emulates  a  large  structure  that  mirrors  cells  [44]  as  shown  in  two  dimensions  in  Figure  3.3. 
Other  periodic  boundary  conditions  emulate  structures  with  other  arrangements  of  re¬ 
peating  cells.  In  the  following  it  is  shown  that  the  boundary  conditions  chosen  have 
negligible  effect  on  energy  localization  and,  by  extrapolation,  on  effective  permittivity. 

3.4.2  Energy  Localization 

EM  remote  probing  of  compounds  presents  the  material  with  an  RF  pulse  resulting 
in  time-localized  energy  concentrations  rather  than  steady-state  energy  concentrations. 
Temporal  localization  of  energy  is  explored  in  this  section  using  the  Gaussian  pulse 
shown  in  Figure  3.4  (with  an  appreciable  frequency  content  from  0  to  50  GHz)  and  the 
irregularly-shaped  crystal  shown  in  Figure  3.1.  An  irregularly-shaped  crystal  of  HMX 


29 


Figure  3.3:  Cross  section  of  the  larger  emulated  structure  based  on  a  unit  cell  shown 
in  Figure  3.1.  The  electric  and  magnetic  boundary  conditions  of  the  TEM  environment 
create  electric  and  magnetic  mirror  boundaries,  extended  here  for  the  25  inner  most  cells. 
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Figure  3.4:  Gaussian  pulse  excitation  signal  used  in  simulation. 

(having  low  real  relative  permittivity  of  2.03  [53],  [54])  in  a  binder  having  a  relative 
permittivity  of  1,  and  in  the  TEM  environment  has  the  (near-maximum)  temporal  energy 
distribution  shown  in  Figure  3.5.  That  is,  this  is  the  EM  energy  distribution  close  to  the 
point  in  the  mixture  that  has  the  highest  electrical  energy  concentration  at  any  position 
and  time. 

It  is  seen  that  the  electric  energy  density  localizes  around  the  edges  and  corners  of  the 
crystal  as  expected  [45].  These  corners  and  edges  are  specifically  the  ones  closest  to  the 
outer  box,  i.e.  the  adjacent  crystal  for  the  emulated  structure,  in  the  polarization  axis 
of  the  applied  electric  field  (the  x  axis  here).  Before  considering  alternative  boundary 
conditions,  consider  energy  localization  for  a  high  permittivity  material.  For  a  zirconia- 
like  high  permittivity  inclusion  (with  real  relative  permittivity  of  28  [55])  the  resulting 
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Figure  3.5:  Peak  electric  energy  density  at  t  =  40  ps  for  the  HMX-based  composite 
(calculated  as  J/m3  and  here  normalized  to  a  maximum  of  1)  with  highest  density  at  the 
edges  and  corners  indicated  by  the  black  oval.  Calculated  using  the  TEM  environment 
with  the  applied  electric  field  polarized  in  the  x  direction  and  propagating  in  the  positive 
z  direction. 
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Figure  3.6:  Peak  electric  energy  density  for  the  zirconia-like  high  permittivity  inclu¬ 
sion  normalized  to  the  energy  density  in  Figure  3.5.  Conditions  correspond  to  those  in 
Figure  3.5. 


electric  energy  density  in  the  composite  is  given  in  Figure  3.6.  The  simulation  conditions 
are  the  same  as  for  Figure  3.5.  Again  the  electric  energy  density  localizes  at  the  same 
corners  and  edges  of  the  inclusion.  Now,  the  peak  electric  energy  density  is  approximately 
8  times  greater  than  with  the  composite  with  the  HMX  (low  permittivity)  inclusion. 

To  investigate  the  effect  of  the  specific  boundary  conditions  on  energy  localization, 
simulations  were  repeated  with  periodic  boundary  conditions  at  the  y-z  and  x-z  planes 
on  the  boundaries.  This  creates  a  repeating  non-mirrored  periodic  structure  in  the  x 
and  y  directions.  With  these  boundary  conditions  it  is  also  necessary  to  change  the 
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Figure  3.7:  Normalized  peak  electric  energy  density  for  the  HMX  (low  permittivity) 
inclusion  with  highest  density  at  the  edges  and  corners  indicated  by  the  black  oval.  Plane 
wave  excitation  and  periodic  boundaries  are  utilized.  The  electric  held  is  polarized  in  the 
x  direction  with  propagation  in  the  positive  z  direction. 


source  excitation  from  a  waveguide  port  to  a  plane  wave.  The  plane  wave  is  linearly 
polarized,  propagating  along  the  z  axis  with  electric  held  polarized  in  the  x  direction. 
The  peak  electric  energy  densities  for  the  low-  and  high-permittivity  inclusions  are  shown 
in  Figure  3.7  and  Figure  3.8  respectively.  In  Figure  3.7  and  Figure  3.8,  the  periodic 
boundary  conditions  result  in  localization  of  electric  energy  density  at  the  corners  and 
edges  of  the  inclusion  and  the  relative  energy  distribution  is  approximately  identical  to 
that  obtained  when  the  mirrored  boundaries  were  used  in  simulation  as  in  Figure  3.5  and 
Figure  3.6. 
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Figure  3.8:  Electric  energy  density  for  the  high  permittivity  inclusion  normalized  to  the 
energy  density  in  Figure  3.7  with  the  same  simulation  conditions  as  in  Figure  3.7. 
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Comparing  the  energy  densities  in  Figure  3.5  and  Figure  3.7,  and  in  Figure  3.6  and 
Figure  3.8,  and  comparing  the  scales  in  Figure  3.6  and  Figure  3.8  enables  the  impor¬ 
tant  conclusion  to  be  reached  that  the  simulation  conditions  have  little  effect  on  energy 
density.  Figure  3.6  and  Figure  3.8  show  the  energy  density  when  the  high  permittivity 
inclusion  is  used  with  the  simulation  reported  in  Figure  3.6  using  mirror  boundaries  and 
Figure  3.8  using  periodic  boundaries.  The  peak  energy  density  differs  by  less  than  5%. 
Thus  the  results  in  Figure  3.5  to  Figure  3.8  show  that  there  is  no  significant  dependence 
of  electric  energy  localization  on  simulation  boundary  conditions.  This  conclusion  was 
reinforced  by  simulations  with  the  crystal  rotated  and  energy  density  again  found  to  have 
little  dependence  on  boundary  conditions.  With  these  results,  the  effective  permittivity 
of  the  compounds  can  be  explored  without  being  concerned  with  the  impact  of  boundary 
condition  choice.  Specifically  TEM  conditions  are  used,  which  enables  the  effective  per¬ 
mittivity  to  be  extracted  using  scattering  parameters  that  can  be  conveniently  extracted 
using  the  TEM  environment. 

3.4.3  Effective  Permittivity  with  HMX  Inclusions 

HMX  has  been  measured  to  have  a  relative  electrical  permittivity  from  1  to  6  GHz  of 
er  =  2.03  —  j0.0035  and  this  is  expected  to  be  constant  up  to  20  GHz  [53],  [54],  The 
flat  permittivity  indicates  that  the  loss  mechanism  is  due  to  dielectric  relaxation  and 
that,  at  least  between  1  and  6  GHz,  there  is  negligible  loss  due  to  material  conductivity. 
The  real  component  of  the  effective  relative  permittivity,  3?(eeff),  at  1  GHz  derived  from 
EM  simulations  for  the  HMX-based  compound  versus  filling  factor  is  compared  to  that 
calculated  using  the  conventional  mixing  theories  in  Figure  3.9. 

For  the  results  in  Figure  3.9  derived  from  EM  simulations,  a  total  of  6  propagation 
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Figure  3.9:  Plots  of  the  real  component  of  the  effective  relative  permittivity  versus  vol¬ 
ume  fraction,  i.e.  filling  factor  q,  for  a  compound  with  irregular  HMX  crystal  inclusions. 
The  squares  identify  the  effective  permittivities  calculated  (using  EM  simulated  results) 
as  the  equally  weighted  average  of  the  effective  permittivities  extracted  for  each  of  the 
three  positive  orthogonal  propagation  directions.  The  plus  signs  identify  the  effective 
permittivities  calculated  similarly  but  now  considering  both  positive  and  negative  propa¬ 
gation  directions.  Also  shown  are  the  effective  permittivities  calculated  using  the  Maxwell 
Garnett  and  Bruggeman  mixing  theories. 
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directions  were  considered  in  turn.  Both  the  forward  and  reverse  propagation  directions 
along  each  of  the  3  coordinate  axes  (x,  y,  and  z),  each  corresponding  to  a  unique  polar¬ 
ization  of  electric  held,  E,  are  used  to  gather  six  independent  results  for  effective  per¬ 
mittivity  according  to  Equation  3.3  through  Equation  3.5.  Figure  3.9  compares  3?(eeff) 
taken  as  the  arithmetic  average  of  either  just  the  positive  propagation  directions,  or  of 
all  six  directions.  In  Figure  3.9,  the  forward  and  reverse  propagation  directions  for  a 
given  polarization  yield  the  same  effective  permittivity  within  1%  and  the  two  averaging 
schemes  yield  results  within  1%.  Thus,  while  individual  permittivity  values  differ  among 
x,  y,  and  z  polarizations,  effective  permittivity  is  independent  of  propagation  direction 
for  the  same  electric  held  polarization.  As  an  additional  test,  the  unit  cell  was  doubled  in 
the  direction  of  propagation,  and  the  calculated  effective  permittivity  was  found  to  differ 
from  the  situation  of  a  single  unit  cell  in  the  propagation  direction  by  less  than  0.05%. 

At  least  for  the  low  permittivity  contrast  material  (as  with  the  HMX  compound  here), 
the  results  shown  in  Figure  3.9  for  the  simulated  effective  permittivity  closely  follow  the 
predictions  of  the  Maxwell  Garnett  and  Bruggeman  mixing  theories  even  for  high  filling 
factors.  At  a  45%  volume  fraction,  q  =  0.45,  the  discrepancy  of  the  four  results  is  less 
than  0.94%.  The  Maxwell  Garnett  and  Bruggeman  mixing  theory  results  were  based  on 
an  assumption  that  the  inclusions  are  spheres  and  that  the  spheres  are  not  close  to  each 
other  (corresponding  to  low  filling  factors  typically  under  10%)  [42],  That  is,  for  a  low 
permittivity  contrast  granular  material,  the  Maxwell  Garnett  and  Bruggeman  mixing 
theories  can  be  used  with  irregular  inclusions  and  with  filling  factors  up  to  45%  within 
1%  error. 
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3.4.4  Effective  Permittivity  with  High  Permittivity  Inclusions 

In  this  section  the  previous  results  for  granular  materials  having  low  permittivity  inclu¬ 
sions  are  compared  to  those  for  a  granular  material  having  a  high  permittivity  contrast 
between  the  inclusions  and  the  matrix.  The  material  examined  has  zirconia-like  crystal 
inclusions.  The  zirconia-like  crystals  have  a  relative  permittivity  of  28,  a  tan  5  (ratio 
of  imaginary  over  real  permittivity)  of  0.0009  (measured  at  1  GHz  [55]),  have  negligible 
conductivity,  and  have  the  same  shapes  as  the  HMX  crystals  to  remove  crystal  shape  as  a 
variable.  Using  the  same  simulation  environment  as  used  previously,  the  calculated  effec¬ 
tive  electrical  permittivities  at  1  GHz  are  compared  to  those  derived  using  the  Maxwell 
Garnett  and  Bruggeman  mixing  theories  in  Figure  3.10.  As  before  the  simulated  effective 
permittivity  is  independent  of  propagation  direction.  Now,  however,  there  is  a  significant 
difference  between  the  effective  permittivities  calculated  and  the  classic  mixing  theories. 
The  discrepancy  is  near  1%  up  to  a  volume  fraction  of  0.05  growing  to  22%  and  61% 
differences  from  the  Maxwell  Garnett  and  Bruggeman  mixing  theories,  respectively,  at 
a  volume  fraction  of  0.45.  This  compares  to  less  than  1%  discrepancy  at  this  volume 
fraction  for  the  lower  permittivity  HMX  compound. 

3.4.5  Summary 

The  Maxwell  Garnett  and  Bruggeman  mixing  theories  are  based  on  spherical  and  regu¬ 
larly  ordered  inclusions  and  are  generally  accepted  as  being  applicable  for  volume  frac¬ 
tions  of  0.1.  The  investigations  here  have  shown  that  with  a  small  permittivity  contrast 
between  the  matrix  material  and  the  inclusion  (here  2.03),  the  Maxwell  Garnett  and 
Bruggeman  mixing  theories  are  accurate  at  high  filling  factors  and  for  irregular  struc¬ 
tures.  However,  with  a  high  permittivity  contrast  between  the  inclusion  and  the  matrix 
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Figure  3.10:  Plots  of  the  real  component  of  the  effective  relative  permittivity  (derived 
from  EM  simulation  results)  versus  filling  factor  for  an  irregular  high  permittivity  inclu¬ 
sion  crystal.  The  squares  indicate  the  equally  weighted  average  of  the  effective  permittiv¬ 
ities  calculated  for  each  of  the  three  positive  orthogonal  propagation  directions.  The  plus 
signs  indicate  the  equally  weighted  average  effective  permittivities  of  the  positive  and 
negative  propagation  directions  (6  total).  The  circles  indicate  the  results  from  using  the 
Frequency  Domain  Solver.  Also  shown  are  the  effective  permittivities  calculated  using 
the  Maxwell  Garnett  and  Bruggeman  mixing  theories. 
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there  is  a  significant  departure  between  the  predictions  of  the  conventional  mixing  the¬ 
ories  and  the  simulated  effective  electrical  permittivity  of  granular  compounds.  From 
Figure  3.10  it  is  also  evident  that  for  irregularly-shaped  inclusions,  high  permittivity 
contrast  between  inclusion  and  matrix  material  produces  significantly  different  results 
from  conventional  mixing  theories  at  filling  factors  above  0.05.  Thus  for  high  filling  fac¬ 
tors  as  seen  in  Figure  3.10,  the  traditional  mixing  theories  do  not  accurately  predict  the 
effective  properties  of  irregularly  shaped  mixtures  when  the  permittivity  contrast  of  the 
inclusion  and  of  the  matrix  is  high.  This  observation  is  compatible  with  previously  re¬ 
ported  research  results  for  high  permittivity  contrast  when  the  inclusions  have  a  regular 
shape  [42],  [44], 

3.5  Modeling  of  Multiple  Crystals 

This  section  continues  the  investigation  of  compounds  with  a  high  permittivity  contrast 
scenario  but  now  using  multiple  cube  crystal  inclusions  in  an  ordered  arrangement  in 
order  to  study  the  properties  of  a  large  volume  with  multiple  inclusions  but  now  with 
a  regular  shape.  The  use  of  a  regular  shape  simplified  the  gridding  used  in  simulation 
and  resulted  in  a  simulation  model  that  fit  in  the  available  160  GB  of  RAM.  An  example 
of  such  a  structure  is  shown  in  Figure  3.11.  In  Figure  3.11,  each  cube  is  0.056  mm  on 
a  side  and  the  side  length  of  the  outer  box  is  1.37  mm  so  that  the  volume  fraction  is 
5%.  For  propagation  in  free  space  at  1  GHz,  a  unit  cell  size  of  1.37  mm  corresponds  to 
approximately  0.0046A  (A  is  the  wavelength  measured  in  mm).  Also  shown  in  Figure  3.11 
are  Ports  1  and  2,  the  waveguide  simulation  ports.  The  side  length  of  each  cube  inclusion 
varies  with  volume  fraction  and  increases  to  approximately  0.15  mm  for  a  0.95  volume 
fraction.  The  inclusions  are  the  high  permittivity  material  and  the  matrix  has  a  relative 
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Figure  3.11:  A  total  of  729  cube  crystals  (inclusions)  arranged  on  a  grid  within  a  free 
space  outer  block.  The  volume  fraction  of  inclusions  here  is  5%. 


permittivity  of  1.  Results  of  the  simulated  effective  permittivity  versus  volume  fraction 
are  given  at  1  GHz  in  Figure  3.12. 

In  Figure  3.12,  the  effective  permittivity  derived  from  simulations  closely  matches  the 
predictions  of  the  Maxwell  Garnett  mixing  theory  for  all  volume  fractions.  Comparing 
Figure  3.12  with  Figure  3.10,  with  the  same  high  permittivity  contrast,  irregular  inclu¬ 
sion  shapes  exhibit  greater  effective  permittivity  differences  from  established  theory  (see 
Figure  3.10)  than  when  symmetric  and  ordered  inclusion  crystals  on  a  grid  are  considered 
(see  Figure  3.12).  It  is  also  possible  that  the  larger  structure,  as  with  the  Figure  3.12 
results,  provides  a  greater  averaging  effect  than  the  thin  layer  of  crystals  considered  with 
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Figure  3.12:  Plots  of  real  effective  relative  permittivity,  3fJ(eeff),  versus  various  filling 
factors  for  high  permittivity  cube  inclusion  crystals  on  a  grid  along  with  various  bounds. 

Figure  3.10.  This  is  examined  further  by  considering  properties  of  a  crystal  while  it  is 
rotated. 


3.6  Rotated  Crystals 

This  section  considers  the  electrical  characteristics  of  crystals  that  are  both  regularly- 
shaped  and  irregularly-shaped  when  the  crystals  are  rotated,  with  permittivity  results 
taken  at  1  GHz.  Effective  permittivity  is  calculated  as  the  inclusion  crystals  are  rotated, 
and  this  behavior  is  compared  with  plots  of  cross-sectional  area.  The  zirconia-like  crys¬ 
tals  having  high  permittivity  contrast  are  considered  since  these  were  found  to  have  the 
greatest  deviation  from  established  theory  previously.  However,  the  crystal  is  now  rotated 
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around  each  of  the  coordinate  axes  (x,  y,  and  z)  to  investigate  the  influence  of  crystal 
orientation  on  the  determination  of  effective  relative  permittivity.  These  are  particularly 
interesting  effects  as  they  form  the  basis  of  a  sensing  modality  with  the  source  of  illumi¬ 
nation  moving  rather  than  rotation  of  the  crystal  inclusions.  Note  that  the  manufacture 
of  many  compound  objects  results  in  aligned  crystals. 

3.6.1  Cube-shaped  Rotated  High  Permittivity  Inclusions 

The  analysis  begins  with  a  simple  rotation  scenario,  a  high  permittivity  cube  crystal 
inclusion  with  a  volume  fraction  of  0.10,  see  Figure  3.13.  In  Figure  3.13,  the  inner  cube  has 
a  side  length  of  0.195  mm  and  the  outer  box  has  a  side  length  of  0.420  mm.  At  1  GHz  a  unit 
cell  size  of  0.420  mm  corresponds  to  approximately  0.0014A  (A  is  the  wavelength  measured 
in  mm)  for  propagation  in  free  space.  The  effect  of  rotation  on  effective  permittivity 
with  the  electric  field  polarized  in  the  x  direction  is  shown  in  Figure  3.14.  Figure  3.14 
indicates  that  when  the  inclusion  rotates  around  the  x  axis  (the  direction  of  electric 
field  polarization),  there  is  negligible  change  in  3?(eeff)  of  the  mixture.  However,  3?(eeff) 
does  vary  with  rotation  of  the  crystal  around  the  y  and  z  axes.  Many  possible  sources 
of  this  rotational  dependency  were  found  with  the  strongest  correlation  being  with  the 
projection  of  the  immediate  cross-sectional  area  of  the  crystals  on  the  plane  transverse 
to  the  direction  of  polarization.  Thus  a  plot  of  cross-sectional  area  fraction  is  shown  in 
Figure  3.15.  The  cross-sectional  area  fraction  with  respect  to  a  given  axis  is  the  ratio 
of  the  cross-sectional  area  of  the  inclusions  projected  in  the  plane  normal  to  the  axis 
to  the  area  of  the  projection  of  the  inclusions  plus  surrounding  material  in  that  plane. 
This  ratio  is  calculated  for  the  first  crystal  layer,  which  is  the  layer  that  interacts  most 
strongly  with  an  illuminating  field.  This  ratio  is  shown  in  Figure  3.15  with  respect  to  the 
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Figure  3.13:  Cube-shaped  high  permittivity  inclusion  crystal  at  10%  volume  fraction  in 
the  TEM  simulation  environment. 

x  axis,  here  the  electric  held  polarization. 

By  comparing  Figure  3.14  with  Figure  3.15,  it  is  evident  that  overall  the  effective 
permittivity  of  the  mixture  with  rotated  cube-shaped  inclusion  follows  the  behavior  of 
the  cross-sectional  area  fraction  graph  when  rotated  about  the  corresponding  axis.  If  the 
cross-sectional  area  fraction  is  constant  (as  is  the  case  for  the  x  axis  rotation  here),  the 
effective  permittivity  is  also  approximately  constant.  However,  from  40°  to  50°  rotation 
around  the  y  and  z  axes,  the  cross-sectional  area  fraction  continues  to  increase  but  the 
effective  permittivity  goes  down.  Outside  this  narrow  range  of  angles  there  is  a  high 
correlation  of  the  cross-sectional  area  and  effective  permittivity. 

Next,  another  cube-shaped  high  permittivity  inclusion  is  used  but  this  time  at  20% 
filling  factor  to  see  how  the  effective  permittivity  versus  rotation  angle  is  altered  when 
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Figure  3.14:  The  real  effective  relative  permittivity,  9fJ(eeff),  of  a  high  permittivity  cube- 
shaped  inclusion  at  a  filling  factor  of  0.10  versus  rotation  angle  around  the  x ,  y,  and  z 
axes.  Here,  the  electric  field  is  polarized  in  the  x  direction. 
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Figure  3.15:  Cross-sectional  area  fraction  from  the  perspective  of  the  x  axis  for  a  10% 
filling  factor.  For  other  filling  factors  the  curves  are  directly  scaled. 
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Figure  3.16:  The  real  effective  relative  permittivity,  3fJ(eeff),  of  a  high  permittivity  cube- 
shaped  inclusion  at  20%  volume  fraction  versus  rotation  angles  around  the  x,  y,  and  z 
axes.  The  electric  field  is  polarized  in  the  x  direction. 

the  filling  factor  is  increased.  The  results  for  effective  permittivity  from  these  simulations 
are  given  in  Figure  3.16.  The  same  behavior  seen  with  lower  filling  factor  is  observed  but 
now  exaggerated. 

For  the  most  part,  effective  permittivity  follows  the  behavior  of  the  cross-sectional 
area  fraction.  Now  however,  there  is  a  wider  region  over  which  the  cross-sectional  area 
fraction  increases  but  the  effective  permittivity  goes  down.  In  Figure  3.16  this  dip  is 
shown  to  be  from  35°  to  55°  rotation  compared  to  40°  to  50°  rotation  for  the  10% 
volume  fraction  simulations  in  Figure  3.14.  For  both  10%  and  20%  volume  fraction,  as 
the  crystal  begins  to  rotate  (y  and  z  axis  rotation  here)  the  cross-sectional  area  fraction 
goes  up.  The  volume  fraction  is  constant  but  because  the  cross-sectional  area  fraction  is 
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going  up  more  of  the  crystal  is  seen  by  the  electric  held  and  the  effective  permittivity 
goes  up.  However,  an  inverse  relationship  develops  between  cross-sectional  area  fraction 
and  effective  permittivity  centered  around  45°.  Since  the  inner  crystal  is  rotating,  at 
45°  rotation  it  is  closest  to  the  nearest  crystal.  The  crystals  getting  close  to  each  other 
causes  another  phenomenon  to  take  over  where  an  inverse  relationship  develops  between 
cross-sectional  area  fraction  and  effective  permittivity. 

The  range  over  which  the  dips  in  Figure  3.14  and  Figure  3.16  occur  can  be  understood 
by  comparing  the  structures  in  these  two  situations.  In  the  20%  volume  fraction  scenario 
in  Figure  3.16,  the  inner  crystal  is  taking  up  more  space  inside  the  outer  box  compared 
to  the  10%  volume  fraction  case  in  Figure  3.14.  In  Figure  3.16  the  region  of  the  dip  has  a 
larger  20°  range  for  the  20%  volume  fraction  scenario  because  (starting  with  0°  rotation 
as  shown  in  Figure  3.13)  it  takes  less  rotation  for  the  crystal  to  get  close  to  the  outer 
box.  There  is  a  smaller  10°  range  of  the  dip  for  the  10%  volume  fraction  scenario  in 
Figure  3.14,  where  because  the  crystal  is  so  small  within  the  outer  box  it  takes  close  to 
45°  rotation  for  the  crystals  to  get  close  enough  to  each  other  for  the  inverse  relationship 
to  develop. 

In  effect,  for  both  the  10%  and  20%  volume  fraction  mixtures  with  cube-shaped 
inclusions  all  of  the  inclusions  are  rotating  and  at  0°  and  90°  the  crystals  are  at  maximum 
separation.  At  45°  the  crystals  are  closest  to  each  other,  where  fringe  effects  influence  the 
effective  permittivity. 

3.6.2  Irregular  Rotated  High  Permittivity  Inclusions 

In  this  section,  the  same  high  permittivity  material  properties  are  used  for  the  inclusion, 
but  this  time  an  irregular  structure  is  used  instead  of  the  cubes  from  the  previous  sec- 
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Figure  3.17:  The  real  effective  relative  permittivity,  9?(eeff),  of  a  high  permittivity 
irregular-shaped  inclusion  at  20%  Filing  factor  versus  rotation  angle  around  the  x,  y, 
and  z  axes.  The  electric  field  is  polarized  in  the  x  direction. 

tion.  From  Sections  3.4  and  3.5,  the  greatest  deviation  from  established  mixing  theories 
occurred  with  high  permittivity  irregular  inclusions,  and  these  structures  are  studied  in 
more  detail  here.  The  shape  of  the  inclusion  is  shown  in  Figure  3.1.  This  allows  a  more 
realistic  and  complicated  structure  to  be  rotated.  For  all  of  the  following  rotations  the 
volume  fraction  is  kept  constant  at  20%.  First,  the  electric  field  is  polarized  in  the  x 
direction,  and  results  of  permittivity  versus  rotation  angle  are  given  in  Figure  3.17. 

In  Figure  3.17,  rotation  around  the  same  axis  as  electric  field  polarization  shows 
virtually  no  deviation  in  effective  permittivity  with  rotation  angle.  However,  rotation 
around  the  other  axes  shows  unique  oscillating  behaviors.  In  an  attempt  to  understand 
the  reasons  for  these  oscillations  in  permittivity,  Figure  3.18  shows  cross-sectional  area 
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Figure  3.18:  Cross-sectional  area  fraction  of  the  rotated  irregular  crystal  from  the  per¬ 
spective  of  the  x  axis  (electric  held  polarization  axis)  as  the  crystal  is  rotated  around  a 
given  axis. 

fraction  looking  along  the  x  axis  (the  electric  held  polarization  axis)  as  the  irregular 
crystal  is  rotated.  By  comparing  Figure  3.17  with  Figure  3.18,  the  behavior  of  cross- 
sectional  area  fraction  with  rotation  follows  that  of  effective  permittivity  for  rotation 
around  the  x  (electric  held  polarization  axis)  and  y  axes.  The  cross  sectional  area  with 
rotation  around  the  z  axis  oscillates  just  as  the  effective  permittivity  does  when  rotated 
around  the  z  axis,  but  with  an  inverse  relationship.  Still,  correlation  between  cross- 
sectional  area  fraction  and  effective  permittivity  is  evident. 

Next,  the  electric  held  is  polarized  in  the  y  direction  and  the  crystal  is  rotated  around 
each  axis.  Results  of  effective  permittivity  are  given  in  Figure  3.19.  Again,  in  Figure  3.19 
there  is  constant  permittivity  when  the  irregular  inclusion  is  rotated  around  the  y  axis 
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Figure  3.19:  The  real  effective  relative  permittivity,  9?(eeff),  of  a  high  permittivity 
irregular-shaped  inclusion  at  20%  filling  factor  versus  rotation  angle  around  the  x,  y, 
and  z  axes.  The  electric  field  is  polarized  in  the  y  direction. 


(electric  field  polarization  axis),  and  oscillating  behavior  with  rotation  around  the  other 
axes.  A  plot  of  the  cross-sectional  area  fraction  is  provided  when  looking  from  the  per¬ 
spective  of  the  y  axis  in  Figure  3.20.  A  comparison  of  Figure  3.19  with  Figure  3.20  shows  a 
direct  correlation  between  cross-sectional  area  fraction  with  rotation  and  effective  permit¬ 
tivity  for  rotation  around  the  x  and  y  (electric  field  polarization)  axes.  Rotation  around 
the  z  axis  causes  cyclical  variation  in  both  figures,  but  with  an  inverse  relationship. 

Analyzing  Figure  3.17  through  Figure  3.20,  it  is  evident  that  there  are  two  situations 
where  the  effective  permittivity  variation  does  not  show  a  direct  relationship  to  the  cross- 
sectional  area  fraction.  These  situations  are  rotation  about  the  z  axis  for  electric  field 
polarized  in  the  x  direction  (Figure  3.17  and  Figure  3.18)  and  rotation  about  the  z  axis 
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Figure  3.20:  Cross-sectional  area  fraction  of  the  rotated  irregular  crystal  from  the  per¬ 
spective  of  the  y  axis  (electric  held  polarization  axis)  as  the  crystal  is  rotated  around  a 
given  axis. 
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for  electric  field  polarized  in  the  y  direction  (Figure  3.19  and  Figure  3.20).  Analyzing 
Figure  3.1,  it  is  seen  that  these  two  situations  involve  the  same  type  of  rotation.  In 
Figure  3.1,  looking  down  the  x  or  y  axes  and  rotating  the  crystal  around  the  z  axis 
causes  the  crystal  to  keep  the  same  width  (in  the  z  direction)  and  rotate  end  over  end. 
In  this  situation  entire  sides  of  the  crystals  are  becoming  close  to  each  other.  The  degree 
of  proximity  is  greater  than  in  the  other  situations.  For  example,  in  rotation  about  the 
y  axis  for  electric  held  polarized  in  the  x  direction  (Figure  3.17),  the  edges  rather  than 
the  faces  of  the  crystals  approach  each  other.  So,  just  as  in  Figure  3.14  and  Figure  3.16 
where  the  cube  crystals  approaching  each  other  resulted  in  an  inverse  relationship  of  the 
effective  permittivity  with  the  cross-sectional  area  fraction,  the  situations  of  rotating  the 
crystal  end  over  end  in  Figure  3.17  through  Figure  3.20  cause  enough  of  the  crystals  to 
be  in  close  proximity  for  an  inverse  relationship  to  develop.  The  simulations  have  shown 
that  close  adjacency  of  neighboring  inclusions  results  in  an  inverse  relationship  between 
cross-sectional  area  fraction  and  effective  permittivity. 

Finally,  the  electric  held  is  polarized  in  the  z  direction  and  the  crystal  is  rotated 
along  each  axis.  The  effective  permittivity  for  this  scenario  is  given  in  Figure  3.21.  Just 
as  shown  in  Figure  3.17  and  Figure  3.19,  when  the  inclusion  is  rotated  around  the  same 
axis  as  the  electric  held  polarization  in  Figure  3.21  (z  axis)  there  is  no  change  in  effective 
permittivity.  Rotation  around  the  x  and  y  axes  shows  the  same  oscillating  behavior.  For 
comparison  purposes,  Figure  3.22  shows  cross-sectional  area  fraction  looking  down  the  z 
axis  (electric  held  polarization  axis)  as  the  irregular  crystal  is  rotated. 

Comparison  between  Figure  3.21  and  Figure  3.22  shows  that  the  cross-sectional  area 
fraction  with  rotation  from  the  perspective  of  the  z  axis  (electric  held  polarization  axis) 
shows  the  same  behavior  as  effective  permittivity  of  the  rotated  inclusion  for  all  axes. 
In  Figure  3.21  and  Figure  3.22,  there  is  no  inverse  relationship  with  cross-sectional  area 
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Figure  3.21:  The  real  effective  relative  permittivity,  5R(eejj),  of  a  high  permittivity 
irregular-shaped  inclusion  at  20%  filling  factor  versus  rotation  angle  around  the  x,  y, 
and  z  axes.  The  electric  field  is  polarized  in  the  z  direction. 
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Figure  3.22:  Cross-sectional  area  fraction  of  the  rotated  irregular  crystal  from  the  per¬ 
spective  of  the  z  axis  (electric  held  polarization  axis)  as  the  crystal  is  rotated  around  a 
given  axis. 
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fraction  because  the  situation  of  rotating  end  over  end  (z  axis  rotation)  with  sides  of  the 
crystals  getting  close  to  each  other  does  not  change  the  cross-sectional  area  fraction  or 
the  effective  permittivity. 

Overall,  rotating  the  inclusion  crystals  provides  strong  evidence  for  a  relationship  be¬ 
tween  the  cross-sectional  area  fraction  of  inclusions  and  effective  permittivity  of  mixtures. 
For  both  regular  and  irregular  crystals,  the  effective  permittivity  and  cross-sectional  area 
fraction  share  a  direct  relationship  except  in  situations  where  a  significant  amount  of  the 
inner  crystal  reaches  close  to  the  outer  box,  or  to  adjacent  crystals,  where  an  inverse 
relationship  develops. 

Rotation  of  the  inclusion  crystals  and  changing  the  electric  held  polarization  axis  high¬ 
lights  the  anisotropy  of  the  composites  studied.  Effective  permittivity  is  found  to  differ 
based  on  electric  held  polarization  direction.  The  results  in  Figure  3.9  and  Figure  3.10 
average  together  the  various  permittivities  obtained  from  varying  the  polarization  di¬ 
rection  to  obtain  a  single  effective  permittivity.  The  sections  describing  rotation  also 
highlight  anisotropy  even  at  a  constant  volume  fraction.  We  change  the  electric  held  po¬ 
larization  axis  and  rotate  the  inclusion  crystal,  and  the  signihcant  differences  observed 
are  highlighted  in  the  figures. 

3.6.3  Impact  of  Operating  Frequency 

The  results  in  the  preceding  sections  present  the  effective  permittivity  at  a  low  frequency 
(1  GHz)  where  the  inclusions  are  small  compared  to  a  wavelength.  To  see  the  impact  on 
effective  permittivity  when  the  operating  frequency  varies,  Figure  3.23  shows  effective 
permittivity  for  the  high  permittivity  irregular  inclusion  at  20%  volume  fraction  with 
electric  held  polarized  in  the  z  direction  (results  from  all  rotations  given  in  Figure  3.21) 
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Figure  3.23:  The  real  effective  relative  permittivity,  3?(eeff),  of  the  high  permittivity 
irregular-shaped  inclusion  at  20%  volume  fraction  (results  for  all  rotation  angles  given  in 
Figure  3.21)  for  0°  rotation  as  a  function  of  frequency. 

but  specifically  for  0°  rotation.  Here  the  unit  cell  has  side  length  of  0.26  mm,  correspond¬ 
ing  to  approximately  0.00087A  at  1  GHz  and  0.043A  at  50  GHz  (A  is  the  wavelength 
measured  in  mm)  for  propagation  in  free  space. 

With  a  crystal  size  of  approximately  0.15  mm,  in  Figure  3.23  the  effective  permittivity 
is  flat  through  25  GHz  before  starting  to  increase  by  less  than  0.5%  up  to  50  GHz,  a  neg¬ 
ligible  change  from  the  static  value.  Similar  results  were  obtained  for  other  polarizations 
and  rotation  amounts  that  also  show  a  flat  response  across  frequency. 
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3.7  Summary 


This  chapter  studied  the  effective  permittivity  and  electric  energy  localization  behavior  of 
crystal-based  compounds  with  high  and  low  contrast  between  the  electrical  permittivity 
of  inclusions  and  the  matrix.  The  results  addressed  the  appropriateness  of  the  simulation 
boundary  conditions,  applicability  of  traditional  mixing  rules,  and  characteristics  that 
can  lead  to  new  sensing  modalities.  In  particular  all  investigations  used  a  temporal  pulse 
corresponding  to  a  remote  probing  scenario. 

For  both  mirrored  and  periodic  boundaries,  electrical  energy  localized  to  the  greatest 
extent  on  the  edges  and  corners  of  the  inclusion  crystal  closest  to  another  crystal  in  the 
direction  of  electric  field  polarization.  The  energy  localization  varied  with  time  indicating 
the  importance  of  not  using  steady-state  analysis  to  determine  peak  energy  density.  In 
particular,  it  is  only  necessary  to  obtain  high  energy  concentration  at  one  point  in  space 
and  time  to  activate  an  energetic  material.  The  energy  distribution  and  temporal  response 
were  virtually  independent  of  the  boundary  conditions  used  in  simulation.  It  is  therefore 
reasonable  to  assume  that  the  boundary  conditions  have  negligible  effect  on  the  effective 
permittivity  of  mixtures  extracted  from  EM  simulations. 

The  effective  permittivities  were  extracted  for  both  forward  and  reverse  propagation 
and  for  each  electric  field  polarization  (x,  y,  and  z).  Since  the  irregular  structures  are  not 
symmetric  with  respect  to  all  polarization  directions,  it  is  not  surprising  that  directional 
dependency  of  the  effective  permittivity  was  observed.  However,  for  the  same  electric 
field  polarization  the  effective  permittivity  did  not  depend  on  forward  or  backward  prop¬ 
agation  direction. 

The  classical  mixing  rules  for  spherical  inclusions,  Maxwell  Garnett  and  Bruggeman, 
were  compared  to  the  effective  permittivities  extracted  from  EM  simulations.  With  low 
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permittivity  contrast  between  inclusion  and  matrix  material  (here  2.03),  there  was  little 
dependence  on  inclusion  shape  for  effective  permittivity  up  to  a  volume  fraction  of  0.45 
(typically  the  highest  volume  fractions  of  energetic  materials  not  of  military  grade).  For 
high  permittivity  contrast  (here  28),  the  effective  permittivity  derived  was  significantly 
different  from  the  prediction  of  the  classic  mixing  theories.  For  ordered  cube  inclusions  on 
a  grid,  the  behavior  of  the  effective  permittivity  derived  from  EM  simulations  was  similar 
to  that  of  the  Maxwell  Garnett  theory  for  all  volume  fractions  (up  to  0.45),  while  neither 
Maxwell  Garnett  nor  Bruggeman  were  able  to  correctly  predict  the  effective  permittivity 
of  mixtures  with  inclusions  having  irregular  shapes  for  volume  fractions  above  0.05. 

One  of  the  purposes  of  the  study  was  to  explore  new  sensing  modalities,  that  is  phys¬ 
ical  behaviors  that  could  be  observed  remotely  with  the  focus  here  being  the  dependence 
of  effective  permittivity  on  observation  direction.  The  effect  of  a  moving  probe  was  emu¬ 
lated  by  rotating  the  crystal  in  the  compound.  When  the  inclusions  were  rotated,  there 
was  a  variation  in  effective  permittivity  of  the  mixture  leading  to  anisotropy  even  at  a 
constant  volume  fraction.  Thus,  the  effective  permittivity  of  a  mixture  varies  as  the  angle 
of  observation  changes.  As  the  crystal  inclusions  were  rotated  (corresponding  to  a  varia¬ 
tion  in  the  angle  of  observation  for  a  mixture  fixed  in  position),  the  effective  permittivity 
for  a  given  electric  field  polarization  axis  was  highly  correlated  to  the  fractional  cross- 
sectional  area  of  the  first  crystal  layer  from  the  perspective  of  the  axis  of  polarization. 
It  is  this  first  crystal  layer  that  most  strongly  interacts  with  an  applied  EM  pulse.  An 
inverse  relationship  between  effective  permittivity  and  cross-sectional  area  fraction  was 
observed.  Note  that  manufactured  energetic  materials  will  typically  have  aligned  crystals. 
That  is,  there  is  a  dependency  of  the  observed  effective  permittivity  on  the  polarization 
and  on  the  angle  of  observation.  This  dependency  is  not  accounted  for  in  the  classic 
Maxwell  Garnett  and  Bruggeman  mixing  theories  and  deserves  further  investigation. 
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The  next  chapter  continues  the  study  of  effective  permittivity  of  composites  and  en¬ 
ergy  localization,  but  increases  the  complexity  of  the  mixtures  that  are  analyzed.  Specifi¬ 
cally,  this  includes  mixtures  with  hundreds  of  disordered  inclusions  both  of  irregular  and 
cube  shapes. 
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Chapter  4 


Electromagnetic  Properties  of 
Mixtures  with  Many  Disordered 
Inclusions 

4.1  Introduction 

The  previous  chapter  analyzed  EM  properties  in  relatively  simple  scenarios  (single  in¬ 
clusions  and  ordered  inclusions  on  a  grid).  To  analyze  more  complex  structures  with 
greater  disorder,  this  chapter  studies  the  EM  properties  of  two-component  mixtures  in¬ 
volving  many  disordered  regularly  and  irregularly  shaped  crystals.  The  effective  relative 
permittivities  are  calculated  utilizing  the  time-domain  Finite  Integration  Technique.  The 
effective  permittivity  of  disordered  mixtures  deviates  from  established  mixing  theories 
especially  in  cases  of  high  permittivity  contrast  between  inclusions  and  matrix  material, 
and  is  strongly  correlated  to  the  cross-sectional  area  of  the  inclusion  crystals.  Electric 
energy  density  localizes  at  the  edges  and  corners  of  inclusions  in  a  manner  independent  of 
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inclusion  shape  and  influenced  by  EM  propagation  direction  and  surrounding  inclusions. 
For  mixtures  with  both  disordered  irregular  and  more  organized  cube  inclusions,  energy 
localization  increases  as  the  EM  signal  travels  through  the  mixture  before  decreasing 
due  to  attenuation  of  the  propagating  EM  signal.  With  a  large  number  of  inclusion  crys¬ 
tals  (here  in  the  hundreds)  it  is  found  that  the  impact  on  effective  permittivity  from 
differences  in  individual  inclusion  shapes  is  negligible. 

Natural  and  manufactured  composites  comprise  two  or  more  components,  and  while 
the  electrical  properties  of  the  individual  components  may  be  well  known,  the  electrical 
properties  of  a  composite  are  generally  thought  to  be  dependent  on  the  shape,  orienta¬ 
tion,  surface  structure,  and  distribution  of  components.  An  understanding  of  the  electrical 
properties  of  component  mixtures,  including  energy  localization,  is  important  in  nonde¬ 
structive  testing,  predicting  the  response  of  rocket  fuel  and  explosives  to  EM  insult, 
remote  sensing,  and  in  industrial  heating  and  curing. 

The  EM  behavior  of  such  a  non-magnetic  complex  object  is  largely  described  by  its 
effective  permittivity.  Studies  of  effective  medium  properties  have  used  various  simulation 
methods.  The  effective  permittivity  of  a  mixture  has  been  calculated  using  finite  difference 
time-domain  EM  simulation  in  two  dimensions  [36],  [45]  and  three  dimensions  [38],  [39], 
in  two  dimensions  using  frequency-domain  finite  element  EM  analysis  [47],  and  in  three 
dimensions  using  the  frequency-domain  finite  difference  [37],  [40],  [41]  and  finite  element 
[42],  [43],  [44]  methods.  The  random  combination  of  variously  shaped  inclusions  with  a 
permittivity  contrast  to  the  embedding  matrix  can  result  in  macroscopic  anisotropy.  The 
standard  electrical  characterization  procedure  is  then  to  average  the  effective  permittivity 
calculated  in  each  of  three  orthogonal  directions  to  obtain  an  overall  effective  permittivity 
of  the  mixture  [42],  [43],  [44],  [56]. 

Earlier  studies  considered  a  single  type  of  inclusion  in  a  matrix  material  with  permit- 
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tivity  contrasts  between  inclusion  and  matrix  of  up  to  10  to  1  [38],  [45].  These  studies 
were  limited  to  a  few  inclusions  because  of  computational  complexity,  and  found  that 
the  effective  permittivity  depended  on  the  shape  of  the  inclusion.  The  previous  chapter 
studied  the  effective  permittivities  of  mixtures  modeled  with  a  single  irregular  crystal 
shape  using  the  FIT  and  found  deviation  from  established  mixing  theories  that  was  espe¬ 
cially  significant  with  high  permittivity  contrast  between  inclusion  and  matrix  material. 
This  chapter  models  mixtures  with  hundreds  of  irregularly  and  regularly  shaped  inclu¬ 
sions  in  a  disordered  arrangement.  The  effective  permittivity  results  are  compared  to  the 
predictions  of  classical  mixing  rules  and  to  each  other.  It  is  shown  here  that  with  many 
inclusions  (in  the  hundreds),  even  when  there  is  relatively  high  permittivity  contast  be¬ 
tween  inclusion  and  matrix  material  of  up  to  28  to  1,  the  impact  of  inclusion  shape  on 
effective  permittivity  is  negligible.  In  addition,  mixtures  with  many  irregular  and  regular 
crystals  are  studied  to  provide  insight  into  how  the  combination  of  individual  crystals 
within  a  larger  mixture  impacts  EM  energy  localization  and  the  creation  of  hotspots. 
Again  it  is  seen  that  inclusion  shape  has  little  effect  on  energy  localization  behavior. 

This  chapter  provides  a  basis  for  understanding  how  pulsed  microwave  signals  can  be 
used  to  characterize  materials  and  for  understanding  how  pulsed  microwave  signals  can 
create  hotspots  in  materials  either  intentionally  or  inadvertently.  For  this  reason,  EM 
analysis  is  performed  in  the  time  domain. 


64 


4.2  Classical  Mixing  Rules 


The  classical  Maxwell  Garnett  mixing  theory  gives  the  effective  permittivity,  eeg,  of  a 
two-component  three-dimensional  mixture  with  spherical  inclusions  as  [21] 


eeff  =  e2  +  3  qe2 


G  —  G 

G  +  2e2  —  <?(g  —  e2) 


(4.1) 


In  Equation  4.1  Ci  is  the  inclusion  permittivity,  e2  is  the  matrix  material  permittivity, 
and  q  is  the  filling  factor  (i.e.,  the  volume  fraction  of  the  inclusions).  The  quasistatic 
assumption  inherent  to  the  development  of  Equation  4.1  is  that  the  inclusion  size  is 
much  smaller  than  the  EM  wavelength  so  that  the  effective  permittivity  is  independent 
of  frequency.  The  Maxwell  Garnett  mixing  theory  also  assumes  spherical  non-touching 
inclusions  that  are  far  apart  [42],  [56]. 

The  classical  Bruggeman  mixing  theory  in  three  dimensions  and  for  spherical  inclu¬ 
sions  is  [22] 


n  G  —  eefj  G  —  ceff 
(!  -  g)_  —  +  9 - - —  =  0. 


(4.2) 


e2  +  2ceff  Ci  +  2eeff 

In  effect,  the  Bruggeman  rule  for  the  effective  permittivity  weights  the  contribution  of 
the  inclusions  by  q  and  that  of  the  matrix  material  by  (1  —  q)  [56],  [57]. 

The  Maxwell  Garnett  and  Bruggeman  mixing  laws  can  be  combined  into  one  uniform 
equation  [37],  [44]: 


Gff  —  G 

Gff  +  2e2  +  zz(eefj  —  e2) 


=  <1 


G  —  e2 

Ci  +  2e2  +  z/(eeff  —  e2) 


(4.3) 


where  v  is  a  parameter  used  to  describe  a  given  mixing  equation.  Maxwell  Garnett  can  be 
obtained  with  v  —  0  and  Bruggeman  with  v  =  2  [37],  [44],  However,  the  classical  mixing 
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rules  are  not  directly  applicable  to  mixtures  with  arbitrarily  shaped  inclusions  [42] ,  [56] . 

The  maximum,  eeff,max,  and  minimum,  eeff,min,  possible  effective  permittivities  of  a 
mixture  are  described  by  the  Wiener  bounds  [36],  [37],  [44]: 


=  qe i  +  (1  -  q)e2 


(4.4) 


^eff,min 


^1^2 


(4.5) 


2  +  (l  -  q)ti 

These  maximum  and  minimum  Wiener  bounds  for  permittivity  correspond  to  capacitors 
in  a  circuit  connected  in  parallel  or  series  respectively  [36],  [37].  These  bounds  are  also 
applicable  when  the  permittivities  of  the  components  are  complex 


4.3  Method  of  Simulation 

Calculated  or  measured  scattering  (S'-)  parameters  can  be  employed  to  find  the  effective 
permittivity  of  a  sample  with  a  finite  thickness  d  [49],  [50],  [51],  [52],  [56].  Here  these 
S-parameters  are  derived  from  a  time-domain  EM  analysis  using  a  Gaussian  excitation 
pulse.  In  particular,  GST  Microwave  Studio  [10],  utilizing  the  time-domain  FIT  with  hex- 
ahedral  meshing,  is  used  for  EM  simulations.  To  confirm  the  accuracy  of  the  simulations, 
the  number  of  mesh  cells  was  increased  until  the  effective  permittivity  converged  to  an 
asymptotic  value.  To  achieve  a  change  in  effective  permittivity  of  less  than  1%,  at  least 
5,000,000  mesh  cells  were  required.  Computation  used  an  80  core  machine  with  160  GB 
of  RAM  and  clocking  at  2.66  GHz. 

The  refractive  index  n  of  a  sample  of  length  d  in  the  propagation  direction  can  be 
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calculated  from  the  S'-parameters  (for  excitation  at  Port  1)  as  [51],  [52],  [56]: 


n  =  ± 


1 

— -arccos 
kd 


1 

2 S^i 


(1  -  Sl  +  si) 


(4.6) 


where  the  free  space  wavenumber  k  =  co/c,  oj  is  the  angular  frequency,  and  c  is  the  speed 
of  light.  The  integer  m  indicates  that  multiple  solutions  are  possible.  Also,  the  relative 
wave  impedance  z  is  defined  as  [51],  [52],  [56]: 


(1  +  S„)2  -  SI 
(1  -SnY-Sl' 


(4.7) 


In  Equation  4.6  and  Equation  4.7  the  S'-parameters  are  normalized  to  the  impedance 
of  free  space,  77.  When  d  and  the  wavelength  A  are  comparable,  obtaining  a  unique 
result  for  n  can  be  difficult.  However  in  this  chapter  d  is  less  than  A/4  at  the  frequency 
analyzed,  1  GHz,  so  a  unique  result  for  n  can  be  found  and  the  default  branch  is  used 
with  m  —  0  and  positive  n  in  Equation  4.6.  With  a  passive  material  is  positive,  so 
in  Equation  4.7  the  positive  branch  is  taken.  These  solutions  for  n  and  z  can  be  used  to 
find  the  unambiguous  effective  permittivity  for  a  nonmagnetic  mixture  as  [51],  [52],  [56]: 


eeff  =  — .  (4.8) 

From  [51],  the  calculations  in  Equation  4.6  and  Equation  4.7  utilize  the  transmission  and 
reflection  coefficients  for  the  mode  that  is  propagating,  which,  for  the  situation  analyzed 
in  this  chapter,  is  TEM  through  50  GHz. 

Before  the  composites  can  be  simulated  in  GST,  their  physical  structure  must  first  be 
defined.  The  next  section  describes  an  automated  process  for  creating  complex  mixtures. 
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4.4  Automated  Creation  of  Complex  Composites 


The  manual  creation  of  complex  three-dimensional  structures  for  use  in  engineering  anal¬ 
ysis  is  a  major  obstacle  to  analyzing  physically  realistic  structures.  A  bias  is  invariably 
imposed  when  a  mixture  is  manually  composed  and  the  structure  rarely  is  representative 
of  the  process  by  which  composites  are  fabricated.  Properties  such  as  packing  density 
and  anisotropies  that  seem  to  easily  occur  in  nature  are  very  difficult  to  obtain  with 
manual  arrangements.  This  section  addresses  the  creation  of  complex  three-dimensional 
mixtures,  comprising  crystals  embedded  in  a  matrix,  for  subsequent  EM  analysis.  The 
physically  realistic  arrangement  of  the  crystals  is  facilitated  by  the  use  of  physics  engine 
software,  specifically  the  Bullet  physics  library,  which  renders  the  realistic  effects  in  ad¬ 
vanced  computer  games.  A  composite  mixture  of  crystals  is  created  by  pouring  a  series 
of  random  crystals  into  a  box  with  the  crystals  bouncing  against  each  other  and  aligning 
just  as  they  do  in  the  real  world.  Higher  packing  densities  are  obtained  than  can  be  rea¬ 
sonably  obtained  with  manual  construction.  The  arrangement  of  the  crystals  obtained 
reflects  the  real  world  alignment  of  asymmetric  crystals.  A  composite  is  created  here  and 
used  with  EM  simulation  software  to  investigate  energy  localization  in  materials. 

Two  of  the  biggest  difficulties  in  the  study  of  random  materials  are  the  inability  of 
humans  to  make  something  truly  random  and  the  manual  writing  of  an  input  language  to 
describe  the  structures.  However,  significant  research  such  as  the  detection  and  analysis 
of  energetic  materials  [1]— [5]  can  be  aided  through  the  automated  creation  and  simulation 
of  complex  crystalline  mixtures.  This  section  uses  the  Bullet  physics  library  [11],  tradi¬ 
tionally  used  in  gaming  software,  to  create  randomized  non-overlapping  irregular  crystal 
shapes  with  various  orientations  and  sizes  that  become  part  of  a  mixture  with  many 
crystals.  This  automated  procedure  can  be  run  many  times  and  can  be  used  to  create 
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complex  3D  crystal  mixtures.  Then  this  structure  is  imported  into  an  EM  simulation  tool 
to  study  how  electric  energy  localizes  in  situations  of  many  irregular  crystals. 

4.4.1  Creating  Multiple  Irregular  Crystal  Structures 

A  two-step  process  is  required  to  perform  EM  simulation  for  mixtures  involving  many  ir¬ 
regular  crystals.  First,  a  structure  must  be  created  (a  nontrivial  process  for  mixtures  with 
many  irregularly  shaped  and  scaled  non-overlapping  crystals)  before  it  can  be  simulated. 

A  model  of  a  composite  structure  is  created  with  physical  simulation  software  using 
the  Bullet  physics  library.  First,  a  set  of  bounding  planes  is  constructed  to  form  an  outer 
box.  Then,  variants  of  different  irregular  crystal  shapes  are  created  at  random  positions 
above  the  box.  An  example  of  two  of  these  differently  shaped  irregular  crystals  is  given 
in  Figure  4.1. 

To  fill  each  outer  box,  approximately  500  individually  scaled  and  rotated  irregular 
crystals  are  created.  For  these  studies  a  total  of  seven  different  crystal  shapes  are  made 
before  they  are  scaled  and  rotated.  Each  crystal  has  a  nominal  size  of  roughly  0.15  mm  on 
each  side,  and  each  dimension  is  scaled  by  a  factor  ranging  from  0.1  to  1  with  a  Gaussian 
distribution,  to  represent  crystals  of  varying  sizes  within  a  mixture.  Even  at  a  scaling 
factor  of  1,  the  size  of  a  crystal  is  much  smaller  than  the  EM  wavelength. 

The  physics  simulation  library  simulates  the  crystals  falling  into  the  box,  rotating, 
pushing  against  each  other  and  being  moved  by  other  crystals  until  they  settle.  The  top  of 
the  box  is  closed  with  another  bounding  plane,  deleting  any  crystals  that  extend  outside 
of  the  box.  At  this  point  none  of  these  irregular  crystals  are  touching.  This  process  of 
utilizing  the  falling  crystals  is  repeated  to  obtain  an  array  of  unique  structures.  Each 
time  the  simulation  is  run  a  different  mixture  is  created.  Figure  4.2  and  Figure  4.3  show 
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Figure  4.1:  An  example  of  two  different  irregularly  shaped  inclusion  crystals,  one  with 
8  sides  (left)  and  one  with  7  sides  (right)  as  created  using  the  Bullet  physics  library. 
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Figure  4.2:  Unique  structure  created  with  550  irregular  crystals  inside  an  outer  box 
(with  19.5%  of  the  outer  box  volume  containing  crystals)  as  created  using  the  Bullet 
physics  library. 

two  of  these  structures.  The  packing  densities  obtained  by  manually  inserting  crystals 
reached  as  high  as  11%,  while  the  packing  densities  from  using  the  Bullet  physics  library 
reached  42%.  This  level  is  representative  of  the  packing  densities  obtained  with  realistic 
structures  [59]. 

After  creating  the  irregular  crystalline  mixtures,  they  can  then  be  imported  into  other 
programs  such  as  CST  Microwave  Studio  for  EM  analysis  using  over  5,000,000  mesh  cells 
and  an  80  core  machine  with  160  GB  of  RAM  and  operating  at  2.66  GHz. 

A  structure  with  many  irregular  crystals  is  created  at  a  certain  volume  fraction  us¬ 
ing  the  Bullet  physics  library  with  output  of  the  physics  library  using  Visual  Basic  for 
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Figure  4.3:  Unique  structure  created  with  500  irregular  crystals  inside  an  outer  box 
(with  21.4%  of  the  outer  box  volume  containing  crystals)  as  created  using  the  Bullet 
physics  library. 
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Figure  4.4:  The  same  mixture  shown  in  Figure  4.2  along  with  waveguide  excitation  ports 
after  it  has  been  imported  into  CST  Microwave  Studio  for  EM  simulation. 

Applications  (VBA)  code.  The  code  is  written  against  CSTs  application  programming 
interface  (API)  then  imported  into  CST  Microwave  Studio  for  EM  analysis.  The  structure 
in  Figure  4.2  is  shown  in  Figure  4.4  but  now  in  the  CST  Microwave  Studio  environment 
along  with  waveguide  excitation  ports  1  and  2  used  in  simulation.  Similarly,  the  complex 
mixture  in  Figure  4.3  as  seen  within  CST  Microwave  Studio  is  given  in  Figure  4.5. 

By  using  CST  Microwave  Studio,  an  EM  pulse  can  be  propagated  through  the  com¬ 
posite  incident  at  either  ports  1  or  2  in  Figure  4.4  and  Figure  4.5.  The  time-domain  solver 
within  CST  is  used  to  propagate  a  Gaussian  pulse  (with  frequency  content  of  0  to  50 
GHz)  through  the  mixtures  in  Figure  4.4  and  Figure  4.5.  Results  of  this  simulation  can 
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Figure  4.5:  The  same  mixture  from  Figure  4.3  with  waveguide  excitation  ports  after  it 
has  been  imported  into  CST  Microwave  Studio  for  EM  simulation. 


74 


be  used  to  study  how  EM  energy  behaves  and  localizes  in  time  when  traveling  through 
these  complex  irregular  mixtures. 

4.4.2  Localization  of  Electromagnetic  Energy 

In  CST,  the  irregular  crystals  are  given  a  relative  permittivity  of  28  surrounded  by  matrix 
material  with  a  relative  permittivity  of  1,  so  that  the  permittivity  contrast  between 
the  crystals  and  the  surrounding  material  is  high.  Following  EM  simulation  the  near- 
maximum  electric  energy  density  obtained  for  the  structure  in  Figure  4.4  is  shown  in 
Figure  4.6,  and  for  the  structure  in  Figure  4.5  in  Figure  4.7.  Both  densities  are  at  the 
specific  time  of  44  ps.  In  Figure  4.6,  the  electric  field  is  polarized  along  the  x  axis  and 
the  EM  wave  propagates  in  the  positive  z  direction  (down).  In  Figure  4.7,  the  electric 
field  is  polarized  along  the  x  axis  and  propagation  is  in  the  negative  z  direction  (up). 

The  combination  of  crystals  in  Figure  4.6  is  different  than  the  combination  of  crystals 
in  Figure  4.7,  since  the  process  of  creating,  scaling,  and  combining  the  crystals  is  per¬ 
formed  independently.  So,  EM  energy  travels  through  a  different  combination  of  crystals 
with  each  EM  simulation  run  and  both  structures  show  high  electric  energy  density.  Lo¬ 
calization  occurring  on  the  edges  and  corners  of  the  crystals  is  present  for  both  Figure  4.6 
and  Figure  4.7,  as  has  been  presented  previously  for  crystals  with  a  fractal  shape  [45]. 

As  detailed  in  this  section,  an  automated  creation  method  of  irregular  crystal  based 
composites  using  the  Bullet  physics  library  has  been  used  to  study  EM  propagation 
through  complex  mixtures.  Electric  energy  was  found  to  localize  on  several  of  these 
irregular  crystals.  While  it  would  be  an  arduous  task  to  create  each  crystal  individually 
by  hand,  the  procedure  described  allows  many  irregular  crystals  to  combine  in  close 
proximity  with  each  other.  Complex  structures  can  therefore  be  created  in  a  reasonable 
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amount  of  time.  Once  the  complex  mixtures  are  formed,  they  can  then  be  imported  into 
other  programs  for  further  analysis. 

4.5  Multiple  Crystal  Mixtures 

Now  that  a  process  to  establish  complex  combinations  of  irregular  crystals  has  been  estab¬ 
lished,  mixtures  involving  many  irregular  crystals  are  created  within  CST  before  their  EM 
properties  are  defined.  In  this  section  the  effective  permittivity  of  various  two-component 
mixtures  is  derived.  The  mixtures  have  variously  shaped  crystal  inclusions  having  Erst 
low  and  then  high  permittivity  contrast  with  the  embedding  matrix.  The  mixtures  in¬ 
volving  many  irregular  and  cube-shaped  crystals  created  in  the  CST  Microwave  Studio 
environment  were  defined  using  an  automatic  procedure  based  on  computer  gaming  soft¬ 
ware  [60].  In  particular,  a  computer  game  was  created  in  which  crystals  (the  inclusions) 
were  poured  into  a  box  and  the  physics  engine  of  the  game  [11]  modeled  the  jostling  and 
bouncing  of  the  crystals  as  they  packed  under  gravity. 

4.5.1  Crystal  Structures  in  CST  Microwave  Studio 

The  irregular  crystals  used  as  inclusions  have  a  maximum  size  of  0.19  mm  as  used  in  [56]. 
A  scaling  factor  in  the  range  of  0.1  to  1  (and  so  the  crystal  sizes  are  much  smaller  than 
a  wavelength  at  50  GHz)  was  then  used  to  randomly  assign  sizes  to  each  crystal.  The 
different  scaling  factors  for  each  crystal  create  a  structure  with  an  array  of  crystal  sizes 
but  each  crystal  has  the  same  shape.  An  example  of  a  complex  mixture  with  irregularly 
shaped  crystal  inclusions  in  the  CST  Microwave  Studio  environment  is  given  in  Figure  4.8. 
Additional  mixtures  were  created  still  exhibiting  a  disordered  arrangement  of  crystals, 
but  this  time  with  crystals  shaped  as  cubes.  An  example  of  such  a  structure  is  shown 
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Figure  4.8:  A  total  of  978  irregularly  shaped  inclusion  crystals  inside  of  a  block  (with 
side  length  of  1.37  mm)  in  a  TEM  simulation  environment.  The  volume  fraction  of  the 
inclusions  is  27.5%. 
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Figure  4.9:  A  total  of  1,144  disordered  cube  inclusion  crystals  (each  with  side  length  of 
0.085  mm)  inside  of  a  vacuum  outer  block  (with  side  length  of  1.37  mm).  The  volume 
fraction  of  inclusions  is  27.5%. 

in  Figure  4.9.  In  both  Figure  4.8  and  Figure  4.9  none  of  the  crystals  are  touching  (a 
separation  that  can  be  controlled)  and  the  crystals  are  surrounded  by  matrix  material 
with  a  relative  permittivity  of  1. 

EM  propagation  through  the  materials  was  modeled  by  placing  the  structures  into  a 
parallel  plate  transverse  EM  (TEM)  environment  [45],  [56].  Waveguide  excitation  Ports 
1  and  2  are  defined  on  the  propagation  axis,  the  z  axis  in  Figure  4.8  and  Figure  4.9, 
and  perfectly  matched  layer  (PML)  boundary  conditions  at  the  x-y  boundaries.  This 
eliminated  reflections  from  the  boundary  planes  back  into  the  TEM  structure.  With 
the  crystal  structure  fixed,  but  with  the  boundary  conditions  changed  appropriately, 
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propagation  in  the  x  and  y  directions  were  also  analyzed.  To  obtain  the  ^-parameters  to 
be  used  in  Equation  4.6  and  Equation  4.7,  measurements  were  de-embedded  to  the  surface 
of  the  cubes  in  Figure  4.8  and  Figure  4.9.  Effective  permittivity  was  then  determined 
using  Equation  4.8  and  a  Gaussian  pulse  excitation  signal  with  frequency  content  up  to 
50  GHz,  as  in  [56]. 

4.5.2  Low  Permittivity  Inclusions 

First,  irregularly  and  cube  shaped  inclusions  of  mica-like  crystals  having  low  relative 
permittivity  were  considered  to  establish  a  lower  permittivity  contrast  scenario  between 
matrix  and  inclusions.  Specifically,  mica  has  a  relative  permittivity  of  5.4  and  a  tan  5  of 
0.0006  measured  at  1  GHz  [55].  Propagation  was  analyzed  for  each  orthogonal  propaga¬ 
tion  axis,  representing  Ex,  Ey,  and  Ez  polarizations.  These  results  were  then  averaged 
to  obtain  a  single  effective  permittivity  (as  described  in  [40],  [41],  [44],  [56]).  Results  of 
effective  permittivity  at  1  GHz  are  compared  to  the  Maxwell  Garnett  and  Bruggeman 
mixing  theories  and  the  Wiener  bounds  in  Figure  4.10. 

Figure  4.10  indicates  that  the  simulated  effective  permittivities  for  the  low  permittiv¬ 
ity  contrast  situation  are  close  to  the  results  predicted  by  Bruggeman  up  to  a  40%  volume 
fraction.  (This  volume  fraction  is  approximately  the  limit  of  what  can  be  obtained  using 
irregularly  shaped  crystals  of  various  sizes  [60].)  The  simulated  effective  permittivities 
for  the  cube  crystals  and  the  irregularly  shaped  crystals  at  the  same  volume  fraction 
differ  by  less  than  0.5%.  This  indicates  that,  with  low  permittivity  contrast,  the  effective 
permittivity  of  mixtures  with  many  inclusions  in  a  disordered  arrangement  has  very  little 
dependence  on  inclusion  shape. 
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Figure  4.10:  Plots  of  real  effective  relative  permittivity,  3fJ(eeff),  versus  various  filling 
factors  for  low  permittivity  irregular  and  cube  inclusion  crystals  along  with  various  mixing 
theories. 
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Figure  4.11:  Plots  of  real  effective  relative  permittivity,  9?(eefr),  as  a  function  of  filling 
factor  for  high  permittivity  irregular  and  cube  inclusion  crystals  along  with  various  mixing 
theories. 


4.5.3  High  Permittivity  Inclusions 

The  above  simulations  were  repeated  with  the  same  structures,  but  this  time  with  in¬ 
clusion  crystals  of  a  higher  permittivity  yielding  a  higher  permittivity  contrast  with  the 
matrix  material.  To  represent  zirconia-like  inclusions  the  crystals  now  have  a  relative 
permittivity  of  28  with  a  tan  5  of  0.0009  taken  at  1  GHz  [55],  [56].  Results  of  effective 
permittivity  at  1  GHz  for  irregularly  shaped  and  cube  inclusions  are  given  in  Figure  4.11. 

In  Figure  4.11  the  simulated  results  fall  within  the  Wiener  bounds  as  expected.  The 
effective  permittivities  for  the  irregular  inclusions  and  cube  inclusions  for  this  high  con¬ 
trast  situation  also  fall  between  the  results  predicted  by  Maxwell  Garnett  and  Bruggeman 
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and  are  close  to  each  other.  Numerically,  of  the  eight  simulated  volume  fractions  for  the 
mixtures  with  cube  inclusions  (as  can  be  seen  in  Figure  4.11),  six  differ  from  the  irregular 
inclusion  results  on  average  by  1%  or  less.  The  greatest  percent  difference  is  2.5%  (with 
an  absolute  difference  in  effective  permittivity  of  approximately  0.07  occuring  at  32.5% 
volume  fraction).  Figure  4.11  shows  that  even  with  a  high  permittivity  contrast  (here 
28  to  1),  the  effective  permittivity  of  composite  mixtures  with  hundreds  of  disordered 
crystals  still  has  very  little  dependence  on  inclusion  shape.  Earlier  research  [38],  [40],  [45] 
has  shown  significant  differences  in  effective  permittivity  between  inclusions  of  different 
shapes  for  single-inclusion  mixtures.  Also,  there  was  no  significant  convergence  in  the 
effective  permittivity  results  when  up  to  eight  of  each  type  of  inclusion  were  analyzed 
[38], 

4.6  Overall  Inclusion  Structure  and  Permittivity 

This  section  explores  in  greater  detail  how  the  overall  arrangement  of  inclusions  influ¬ 
ences  the  individual  extracted  effective  permittivity  components  (i.e.,  eeff  x,  eesiV,  and 
eefr,z  components  in  the  x,  y,  and  z  directions  respectively)  of  a  mixture. 

4.6.1  High  Permittivity  Inclusions 

The  scenarios  that  lead  to  the  maximum  and  minimum  possible  effective  permittivity, 
as  described  by  the  Wiener  bounds,  are  when  the  inclusion  material  has  the  shape  of  a 
plate  [36].  For  example,  a  structure  having  an  effective  permittivity  corresponding  to  the 
upper  Weiner  bound  is  given  in  Figure  4.12. 

In  Figure  4.12,  the  side  length  of  the  outer  box  is  0.91  mm  and  the  thickness  of  the 
plate  in  the  z  direction  is  0.10  mm,  giving  a  volume  fraction  of  11.0%.  The  Wiener  bounds 
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Figure  4.12:  Structure  showing  a  high  permittivity  inclusion  arranged  as  a  plate  in 
vacuum  with  11.0%  volume  fraction. 
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from  Equation  4.4  and  Equation  4.5  can  be  used  to  find  the  maximum  and  minimum 
effective  permittivity  possible.  For  the  structure  in  Figure  4.12,  these  bounds  are  found 
to  be  3.97  for  the  maximum  9ft(eeff)  and  1.12,  changing  the  propagation  direction,  for  the 
minimum  9ft  (eeg). 

These  results  will  now  be  compared  to  the  simulated  effective  permittivities.  For  the 
structure  in  Figure  4.12  the  extracted  effective  permittivity  depends  on  the  axis  of  the 
electric  field  polarization.  The  structure  in  Figure  4.12  has  the  same  electrical  properties 
for  x-  and  ^/-directed  electric  fields.  For  x  and  y  polarization  3ft(eefr)  from  simulation  was 
found  to  be  identical,  i.e.  3.97.  Notably,  this  corresponds  to  the  maximum  Wiener  bound. 
For  ^-directed  electric  field,  9ft(eeff)  from  simulation  was  found  to  be  1.13,  nearly  identical 
to  the  minimum  Wiener  bound. 

So,  as  described  in  [36],  these  results  confirm  that  the  maximum  effective  permittivity 
at  a  given  volume  fraction  occurs  when  the  largest  face  of  the  inclusion  plate  is  parallel  to 
the  electrical  polarization  of  the  TEM  field.  The  minimum  effective  permittivity  occurs 
when  the  largest  face  of  the  inclusion  plate  is  perpendicular  to  the  electrical  polarization 
of  the  TEM  field. 

4.6.2  Relationship  to  Cross-Sectional  Area 

In  an  effort  to  relate  the  arrangement  of  the  plate  in  Figure  4.12  to  the  maximum  or  min¬ 
imum  effective  permittivity,  the  cross-sectional  areas  of  the  inclusion  from  Figure  4.12 
are  studied.  Analyzing  Figure  4.12,  the  maximum  effective  permittivity  occurs  when  the 
inclusion  has  the  minimum  possible  cross-sectional  area  normal  to  the  direction  of  the 
electric  field  polarization.  In  this  scenario,  the  narrow  edge  of  the  plate  is  ‘seen’  first  by  the 
electric  field.  Conversely,  the  minimum  effective  permittivity  occurs  when  the  inclusion 
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has  the  maximum  possible  cross-sectional  area  normal  to  the  electric  field  polarization 
direction.  Here,  the  wide  face  of  the  plate  is  presented  to  the  incident  electric  field.  So, 
this  suggests  that,  for  an  inclusion  shaped  as  a  plate,  there  is  an  inverse  relationship 
between  the  presented  cross-sectional  area  of  the  inclusion  and  the  effective  permittivity 
of  the  mixture.  This  conclusion  is  consistent  with  [56]  where  a  similar  inverse  relation¬ 
ship  between  effective  permittivity  and  cross-sectional  area  was  developed  for  a  single 
arbitrarily  shaped  inclusion  crystal. 

To  investigate  the  cross-sectional  area  relationship  for  structures  containing  many 
inclusions,  a  mixture  containing  irregular  crystals  with  a  15%  volume  fraction,  as  shown 
in  Figure  4.13,  and  a  mixture  containing  cube  crystals  also  with  a  15%  volume  fraction, 
as  shown  in  Figure  4.14,  were  considered.  Cross-sectional  area  was  calculated  using  rays 
cast  along  the  appropriate  axis  using  the  Bullet  physics  library  [11],  For  example,  to 
determine  the  cross-sectional  area  presented  by  crystals  from  the  perspective  of  the  x 
axis,  the  y-z  face  of  the  outer  box  is  divided  into  an  N  x  N  grid  (here  N  =  500).  For 
every  point  on  the  grid,  a  ray  is  cast  perpendicularly  to  the  y-z  face  (in  the  x  direction) 
from  one  side  of  the  outer  box  to  the  other.  If  the  ray  hits  a  crystal,  a  hit  is  recorded. 
The  cross-sectional  area  of  the  crystals,  %Crystais,  is  then  defined  as: 

^crystals  =  M  Acen,  (4.9) 

where  M  is  the  number  of  hits  and  %cen  is  the  area  of  one  grid  cell.  The  cross-sectional 
area  fraction,  AWl  (w=x,  y  or  z)  is  further  defined  for  a  given  axis  perspective  (the  x 
direction  is  used  as  an  example  here)  as: 

Ax  =  (4.10) 
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Figure  4.13:  Structure  showing  894  high  permittivity  irregular  inclusions  inside  a  vac¬ 
uum  outer  block  (with  side  length  1.46  mm).  The  volume  fraction  of  inclusions  is  15%. 


Figure  4.14:  Structure  showing  624  high  permittivity  cube  inclusions  (each  with  side 
length  0.085  mm)  inside  a  vacuum  outer  block  (with  side  length  1.37  mm).  The  volume 
fraction  of  inclusions  is  15%. 
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Figure  4.15:  Real  effective  permittivity,  9fJ(eeff),  plotted  as  a  function  of  cross-sectional 
area  fraction,  Aw.  The  volume  fraction  of  inclusions  is  15%.  The  letters  near  the  data 
points  indicate  the  electric  field  polarization  axis  for  each  type  of  crystal  as  well  as  the 
cross-sectional  area  fraction  direction. 


where  %tota l  is  the  total  area  of  the  outer  box.  The  cross-sectional  area  of  crystals  and 
cross-sectional  area  fraction  are  calculated  in  this  manner  for  the  x ,  y,  and  z  directions 
independently. 

With  high  permittivity  inclusions  at  15%  volume  fraction,  the  effective  permittivi¬ 
ties  and  corresponding  cross-sectional  area  fractions,  Aw,  for  each  axis  are  given  in  Fig¬ 
ure  4.15.  Comparing  all  cases  presented  in  Figure  4.15,  the  highest  cross-sectional  area 
fraction  value  corresponds  to  the  lowest  effective  permittivity  component  and  vice  versa, 
suggesting  an  inverse  relationship. 

The  differences  in  magnitude  of  the  effective  permittivity  among  the  polarizations  as 
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presented  in  Figure  4.15  can  be  qualitatively  related  to  the  cross-sectional  area  fraction 
differences  by  analyzing  the  mixtures  in  Figure  4.13  and  Figure  4.14.  Based  on  the  way 
the  irregular  structures  are  created  using  the  Bullet  physics  library,  the  crystals  fall  from 
above  (with  respect  to  the  y  axis)  and  fall  to  the  bottom  of  the  outer  box  before  settling. 
The  crystals  therefore  become  tightly  packed  on  the  bottom  of  the  box  as  can  be  seen  in 
Figure  4.13  and  Figure  4.14. 

This  means  that  from  the  perspective  of  the  y  axis,  the  inclusion  crystals  have  a  high 
cross-sectional  area  fraction,  Ay.  From  the  perspective  of  the  x  and  z  axes,  the  cross- 
sectional  areas  of  the  inclusions  are  relatively  small,  and  so  Ax  and  Az  are  smaller.  This 
can  also  be  seen  in  Figure  4.13  and  Figure  4.14,  where  from  the  perspective  of  the  x  and 
z  axes  the  upper  part  of  the  box  is  the  background  material  (with  relative  permittivity 
of  1). 

Relating  this  to  the  permittivity  values,  a  high  Aw  from  the  perspective  of  the  TEM 
electric  held  polarization  axis,  w,  leads  to  a  low  effective  permittivity  (e.g.  here  Ay  is  high 
and  the  effective  relative  permittivity  component  in  the  y  direction,  eeff)3/,  is  low)  while 
the  other  polarizations  have  lower  Aw  and  higher  effective  permittivity.  This  is  the  same 
inverse  relationship  demonstrated  by  the  plate  of  Figure  4.12.  Since  the  crystals  do  not 
form  a  complete  plate  in  any  direction,  the  maximum  or  minimum  possible  permittivity 
values  as  given  by  the  Wiener  bounds  are  not  reached.  However,  based  on  the  way  the 
crystals  are  packed  in  the  box  in  Figure  4.13  and  Figure  4.14,  they  are  quite  dense 
and  close  to  the  sides  of  the  box.  In  this  situation  it  can  be  said  that  the  inclusions 
exhibit  plate-like  behavior.  This  can  be  seen  through  the  sharing  of  an  inverse  relationship 
between  cross-sectional  area  fraction  and  effective  permittivity  for  both  the  plate  and  the 
densely  packed  crystals.  This  correspondence  is  seen  for  all  of  the  simulation  results  for 
low  and  high  permittivity  inclusions  presented  in  Figure  4.10  and  Figure  4.11  respectively. 
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Since  all  of  the  simulations  covering  the  various  volume  fractions  in  Figure  4.10  and 
Figure  4.11  were  created  in  the  same  manner  (the  crystals  falling  into  the  box  along  the 
y  axis  and  then  settling)  Aw  for  all  simulations  was  greatest  looking  along  the  y  axis  (i.e. 
w=y),  and  Ey  polarization  always  yielded  the  smallest  effective  permittivity  component. 

4.7  Energy  Localization 

The  other  important  phenomenon  that  characterizes  the  electrical  response  to  a  pulsed 
microwave  signal  is  energy  localization,  in  particular  the  concentration  of  electric  energy 
in  a  nonmagnetic  composite.  Composites  respond  to  temporal  energy  localization  by 
delaminating  and  changing  chemical  phase.  Also,  identifying  energy  localization  can  be 
used  as  a  diagnostic  tool.  Many  of  these  effects  respond  over  a  short  time  interval  so  that 
it  is  the  transient  time-domain  EM  response  that  is  more  important  than  the  steady-state 
response  provided  by  a  frequency- domain  EM  analysis.  Energy  has  been  found  to  localize 
on  the  edges  and  corners  of  single  inclusions  [45],  [56].  This  section  explores  localization 
phenomena  taken  using  time-domain  monitors  for  mixtures  containing  many  crystals 
close  together,  highlighting  phenomena  neglected  by  steady-state  analyses. 

4.7.1  High  Permittivity  Inclusions 

First,  electric  energy  density  is  plotted  for  a  complex  mixture  with  32.5%  volume  fraction. 
An  example  showing  near-maximum  electric  energy  density  in  a  high  permittivity  (with 
a  real  relative  permittivity  of  28)  inclusion  mixture  for  electric  field  polarized  along  the 
x  axis  and  positive  z  propagation  direction  is  given  in  Figure  4.16.  The  excitation  signal 
is  a  Gaussian  pulse  with  frequency  content  up  to  50  GHz.  The  side  length  of  the  outer 
box  in  Figure  4.16  is  1.37  mm.  The  bottom  right  black  circle  shown  in  Figure  4.16  has  a 
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Figure  4.16:  Near-maximum  electric  energy  density  showing  localization  at  inclusion 
corners  and  edges,  with  highest  localization  in  the  bottom  two  circles.  The  electric  held 
is  polarized  along  the  x  axis  and  the  propagation  direction  is  positive  z  (down).  The 
volume  fraction  is  32.5%. 


diameter  of  approximately  0.12  mm.  For  reference,  in  vacuum  the  wavelength  at  1  GHz 
is  300  mm  and  at  50  GHz  is  6  mm.  With  er  =  28,  the  wavelength  at  1  GHz  is  56.69  mm 
and  at  50  GHz  is  1.13  mm.  From  Figure  4.11,  the  effective  permittivity  for  this  composite 
was  calculated  as  2.97.  With  a  relative  permittivity  of  2.97,  the  wavelength  at  1  GHz  is 
174.08  mm  and  at  50  GHz  is  3.48  mm.  Comparing  these  sizes,  the  energy  localization 
shown  in  Figure  4.16  occurs  on  the  order  of  0.01A  or  less,  a  fine  resolution  in  terms  of 
wavelength. 

The  near-maximum  energy  density  is  shown  in  Figure  4.16  at  48.1  ps.  The  energy 
localizes  to  the  greatest  extent  on  the  edges  and  corners  of  the  irregular  inclusions  shown 
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Figure  4.17:  Electric  energy  density  with  the  same  conditions  as  in  Figure  4.16,  but  with 
propagation  in  the  reverse  direction  (up).  The  strong  localization  from  Figure  4.16  is  not 
seen  in  the  bottom  two  circles.  Instead,  there  is  higher  localization  in  the  top  circle. 


inside  the  bottom  two  circles  in  Figure  4.16  after  the  signal  has  passed  through  many 
other  inclusions.  Propagation  in  the  reverse  direction,  along  the  negative  z  axis  (up), 
results  in  the  near-maximum  electric  energy  density  scenario  in  Figure  4.17.  This  image 
is  taken  at  the  same  time  and  plane  as  the  structure  shown  in  Figure  4.16.  The  only 
difference  in  simulation  conditions  between  Figure  4.16  and  Figure  4.17  is  an  opposite 
direction  of  propagation. 

The  crystal  indicated  by  the  upper  circle  in  Figure  4.17  shows  greater  energy  density 
compared  to  the  same  location  in  Figure  4.16.  Conversely,  the  energy  density  on  the 
crystals  indicated  by  the  bottom  two  circles  is  smaller.  Comparing  Figure  4.16  and  Fig¬ 
ure  4.17,  the  electric  energy  density  localizes  with  different  magnitudes  depending  on  the 
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direction  of  propagation  of  the  EM  signal.  At  specific  locations,  indicated  by  the  black 
circles  in  Figure  4.16  and  Figure  4.17,  the  magnitude  of  the  energy  density  at  the  edges 
and  corners  is  higher  when  the  signal  has  traveled  through  more  scatterers. 

As  an  additional  test,  the  Gaussian  pulse  duration  was  extended,  with  the  frequency 
content  ranging  from  0  to  25  GHz.  The  energy  was  seen  to  localize  in  the  same  locations 
with  the  same  localization  behavior  as  for  the  0  to  50  GHz  bandwidth  excitation  used 
in  Figure  4.16  and  Figure  4.17.  Halving  the  frequency  range  allowed  the  Gaussian  pulse 
duration  to  double  from  71  ps  to  142  ps.  Comparing  the  peak  localization  magnitudes 
with  Port  2  excitation,  for  the  0  to  50  GHz  excitation  the  maximum  monitored  energy 
density  was  0.0768  J/nr3,  while  for  the  longer  0  to  25  GHz  excitation  the  maximum  was 
0.0992  J/m3.  Overall,  this  result  suggests  that  a  longer  excitation  time  can  result  in  a 
higher  peak  localization  magnitude. 

For  the  situations  represented  in  Figure  4.16  and  Figure  4.17,  the  peak  input  power  is 
1  W.  This  is  input  over  the  entire  port  face,  which  has  a  side  length  of  1.37  mm.  Assuming 
this  1  W  is  input  for  the  0  to  50  GHz  Gaussian  pulse  duration,  71  ps,  means  that  the  total 
input  power  is  7.1  xl0~n  J.  In  order  to  find  an  energy  density  entering  the  structure,  it  is 
necessary  to  find  the  wavelength  using  the  calculated  effective  permittivity  of  2.97  so  that 
the  effective  group  velocity  is  1.74xl08  m/s.  Since  the  Gaussian  pulse  duration  is  71  ps, 
the  length  of  the  pulse  in  the  medium  is  0.012  m.  Multiplying  this  by  the  input  port  area 
gives  an  input  volume  of  2.25 xlO-8  nr3.  Then,  the  energy  density  entering  is  the  input 
power  divided  by  this  volume,  i.e.  0.0032  J/m3.  Comparing  this  value  to  the  peak  energy 
density  of  0.0768  J/m3  in  Figure  4.16,  this  means  that  here  there  is  an  energy  density 
increase  of  24  times  represented  by  this  energy  localization. 

The  physical  interpretation  is  that  the  energy  localization  maximum  occurs  when 
signals  following  multiple  scattering  paths  coalesce  at  the  same  time  and  position.  Such 
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a  hotspot  is  unlikely  to  occur  at  the  surface  of  a  composite  but  within  the  material 
after  multiple  scattering  events  have  occurred  but  before  the  signal  has  been  significantly 
attenuated. 

4.7.2  Electric  Fields  for  Irregular  Inclusions 

To  gain  a  better  understanding  of  the  localization  behavior  occurring  in  the  circles  in 
Figure  4.16  and  Figure  4.17,  simulation  probes  were  inserted  to  measure  the  electric  field 
magnitudes  over  the  course  of  the  simulation.  The  electric  fields  at  the  location  indicated 
by  the  bottom  right  circles  in  Figure  4.16  and  Figure  4.17  are  given  in  Figure  4.18,  at  the 
bottom  left  circles  are  given  in  Figure  4.19,  and  at  the  top  circles  are  given  in  Figure  4.20. 
In  Figure  4.18  through  Figure  4.20,  the  electric  field  localizes  at  specific  crystal  corners 
and  edges  to  a  greater  amount  when  the  EM  signal  has  previously  traveled  through  more 
scatterers.  From  Figure  4.18  and  Figure  4.19,  peak  energy  localization  at  the  bottom 
crystals  occurs  after  the  signal  has  traveled  from  top  to  bottom.  The  energy  density 
for  this  positive  z  propagation  direction  is  shown  in  Figure  4.16  at  the  time  indicated 
in  Figure  4.18  and  Figure  4.19,  i.e.  48.1  ps.  From  Figure  4.20,  peak  localization  at  the 
top  crystal  occurs  after  the  signal  travels  from  bottom  to  top,  with  energy  density  for 
this  negative  z  propagation  direction  at  the  time  indicated  in  Figure  4.20,  48.1  ps,  given 
in  Figure  4.17.  The  electric  field  probes  from  Figure  4.18  through  Figure  4.20  confirm 
the  results  of  energy  density  from  Figure  4.16  and  Figure  4.17  that  traveling  through  a 
greater  number  of  inclusions  increases  localization  magnitude  at  a  specific  location.  Next, 
this  behavior  is  studied  for  a  simpler  case  of  identical  inclusions  placed  on  an  organized 
grid. 
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Figure  4.18:  Comparison  of  the  electric  field  at  the  crystal  corner  indicated  by  the  bot¬ 
tom  right  circle  in  Figure  4.16  and  Figure  4.17  with  propagation  direction.  At  48.1  ps 
there  is  a  significant  difference  in  the  electric  field  and  electric  energy  density  depending 
on  propagation  direction.  Higher  peak  localization  is  seen  with  the  positive  propaga¬ 
tion  direction  (down)  when  the  EM  signal  has  traveled  through  more  of  the  irregular 
surrounding  material. 
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Figure  4.19:  Comparison  of  the  electric  field  at  the  crystal  corner  indicated  by  the  bot¬ 
tom  left  circle  in  Figure  4.16  and  Figure  4.17  with  propagation  direction.  Higher  peak 
localization  is  seen  at  48.1  ps  with  the  positive  propagation  direction  (down)  when  the 
EM  signal  has  traveled  through  more  of  the  irregular  surrounding  material. 
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Figure  4.20:  Comparison  of  electric  field  for  the  top  circle  in  Figure  4.16  and  Figure  4.17 
with  propagation  direction.  At  48.1  ps  there  is  a  significant  difference  in  the  electric  field 
and  electric  energy  density  depending  on  propagation  direction.  Higher  peak  localization 
is  seen  with  the  negative  propagation  direction  (up)  when  the  EM  signal  has  traveled 
through  more  of  the  irregular  surrounding  material. 
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Figure  4.21:  Cube-shaped  inclusions  arranged  in  two  planes.  On  the  right  are  the  isolated 
cubes  and  on  the  left  are  the  back  cubes. 

4.7.3  Behavior  for  Cube-shaped  Inclusions 

From  Section  4.7.2,  the  electric  held  localization  magnitude  was  seen  to  increase  when 
the  EM  signal  traveled  through  a  greater  number  of  scatterers.  This  section  studies  elec¬ 
tric  held  localization  for  a  less  complicated  mixture  with  high  permittivity  cube-shaped 
inclusions  on  a  regular  three-dimensional  grid.  The  aim  of  this  study  is  to  determine  if  a 
simpler  simulation  environment  yields  similar  energy  localization  results.  First,  two  par¬ 
allel  layers  of  cubes  in  the  x-y  plane  are  created  as  shown  in  Figure  4.21.  Each  cube  has 
a  side  length  of  0.1  mm.  The  left  x-y  plane  in  Figure  4.21  has  20  cubes  placed  randomly 
within  the  plane,  two  of  which  are  referred  to  as  the  back  cubes.  After  the  cubes  have 
been  placed  in  the  plane,  the  back  cubes  (the  two  closest  cubes  in  this  plane)  are  then 
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Figure  4.22:  Layers  of  cube-shaped  inclusions  are  added  between  the  layers  in  Figure  4.21 
with  an  increasing  number  of  inclusions  in  each  layer. 

shifted  in  the  positive  z  direction  with  the  same  x  and  y  coordinates  in  order  to  maintain 
the  same  distance  and  orientation  between  the  cubes  but  to  isolate  them  in  their  own 
x-y  plane.  This  pair  of  cubes,  to  be  referred  to  as  the  isolated  cubes,  are  shown  on  the 
right  side  of  Figure  4.21.  The  two  cubes  are  isolated  in  order  to  study  the  influence  of 
surrounding  cubes  within  a  plane  on  energy  localization. 

Next,  additional  cubes  are  placed  in  planes  evenly  spaced  between  those  shown  in 
Figure  4.21.  The  resulting  structure  is  shown  in  Figure  4.22.  The  distance  between  cubes 
in  two  adjacent  x-y  planes  in  Figure  4.22  is  0.05  mm  and  none  of  the  cubes  are  touching. 
Behind  the  two  isolated  cubes  on  the  right,  going  in  the  negative  z  direction  (left),  the 
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Figure  4.23:  The  structure  in  Figure  4.22  (surrounded  by  vacuum)  along  with  waveguide 
ports  for  simulation. 


next  x-y  plane  has  1  cube,  then  3,  5,  10,  and  15  before  the  left-most  plane  with  20  cubes 
is  reached.  Adding  cubes  in  this  manner  gradually  increases  the  number  of  inclusions 
between  the  layers  shown  in  Figure  4.21.  The  goal  is  to  confirm  the  hypothesis  that  EM 
propagation  through  more  inclusions  leads  to  an  increase  in  the  peak  hotspot  magnitude. 

The  total  structure  with  all  cubes  along  with  the  waveguide  Ports  1  and  2  is  given 
in  Figure  4.23.  In  Figure  4.23,  the  isolated  cubes  from  Figure  4.21  and  Figure  4.22  are 
closest  to  Port  1  and  the  back  cubes  are  closest  to  Port  2.  The  outer  box  in  Figure  4.23 
spans  0.8  mm  in  the  x  and  y  directions.  The  two  ports  are  the  same  distance  to  the 
nearest  high  permittivity  cube  in  either  direction,  specifically  1.55  mm.  The  electric  field 
is  polarized  along  the  x  axis  and  z  is  the  propagation  axis.  The  excitation  signal  is  a  0 
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Figure  4.24:  Near-peak  electric  energy  density  for  the  back  layer  of  cubes  at  44  ps  (time 
of  near-maximum  energy  density)  with  highest  magnitude  at  the  edges  and  corners  of 
the  back  cubes  as  indicated  by  the  black  circle.  This  occurs  for  Port  1  excitation. 


to  50  GHz  Gaussian  pulse.  The  near-peak  energy  density  location  in  the  structure  from 
Figure  4.23  is  shown  in  Figure  4.24. 

In  Figure  4.24,  the  highest  energy  density  among  all  the  cubes  occurs  at  the  edges 
and  corners  of  the  back  cubes  with  Port  1  excitation.  Again,  localization  magnitude  is 
greater  after  the  signal  has  traveled  from  Port  1  to  Port  2  through  a  greater  number  of 
scatterers.  The  plane  of  cubes  shown  in  Figure  4.24  is  the  back  layer  with  the  highest 
number  of  cubes  (20)  closest  to  Port  2.  The  back  cubes  indicated  by  the  black  circle  in 
Figure  4.24  have  edges  and  corners  closest  to  each  other  within  that  plane,  and  no  edges 
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Figure  4.25:  Electric  field  magnitudes  for  four  different  locations  in  the  structure  in 
Figure  4.23  as  a  function  of  time.  Peak  localization  magnitude  is  reached  for  the  back 
cubes  with  Port  1  excitation  after  the  signal  has  traveled  through  other  scatterers. 


between  other  adjacent  cubes  are  closer  in  the  entire  structure.  Overall,  adjacent  cubes 
coming  closer  together  also  increases  temporal  energy  localization. 

Electric  field  probes  can  be  used  to  gain  an  understanding  of  how  the  field  magnitude 
changes  over  time.  Probes  are  utilized  to  measure  the  electric  field  at  the  edge  of  one  of 
the  back  cubes  (encompassed  by  the  circle  in  Figure  4.24)  as  well  as  the  same  relative 
location  but  on  the  isolated  pair  of  cubes  closest  to  Port  1  (those  indicated  in  Figure  4.23). 
The  results  from  these  two  probes,  each  probe  measuring  electric  field  with  Port  1  and 
Port  2  excitation  independently,  are  summarized  in  Figure  4.25.  The  isolated  cubes  and 
back  cubes  referenced  in  Figure  4.25  are  analyzed  because  they  represent  the  maximum 
EM  localization  within  the  mixture  shown  in  Figure  4.23.  In  Figure  4.25,  electric  field 
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localization  magnitudes  are  plotted  as  a  function  of  time.  As  shown  in  Figure  4.21,  there 
are  two  pairs  of  cubes  with  the  same  orientation  to  each  other:  the  isolated  cubes  and  the 
back  cubes.  One  pair  is  surrounded  by  many  other  inclusions  in  the  same  plane  (the  back 
cubes  shown  in  Figure  4.21,  in  the  plane  of  cubes  closest  to  Port  2  in  Figure  4.23,  and 
referenced  in  Figure  4.25)  and  the  other  pair  is  isolated  in  a  plane  by  itself  (the  isolated 
cubes  shown  in  Figure  4.21,  in  the  plane  of  cubes  closest  to  Port  1  in  Figure  4.23,  and 
referenced  in  Figure  4.25)  to  study  the  impact  of  surrounding  inclusions.  The  two  cubes 
within  a  pair  are  the  same  distance  from  each  other  and  differ  only  in  their  z  coordinates. 

Analyzing  Figure  4.25,  for  the  back  cubes  closely  surrounded  by  other  high  permittiv¬ 
ity  cubes,  the  peak  electric  held  magnitude  is  higher  for  both  propagation  directions  (i.e. 
for  Port  1  and  Port  2  excitation,  representing  propagation  in  the  negative  and  positive  z 
directions  respectively)  than  for  either  propagation  direction  for  the  isolated  cubes  in  the 
front.  There  are  more  scatterers  in  the  back  layer,  so  the  localization  magnitude  is  higher 
than  in  the  front  where  in  that  layer  there  are  few  scatterers.  In  summary,  a  greater 
number  of  scatterers  within  a  layer  increases  energy  localization. 

For  the  back  cubes  in  Figure  4.25  (in  the  layer  closest  to  Port  2  in  Figure  4.23  and 
Figure  4.24),  looking  at  the  differences  caused  by  changing  the  propagation  direction, 
there  is  higher  localization  for  Port  1  excitation  than  for  Port  2  excitation.  With  excitation 
at  Port  1,  the  signal  goes  through  six  layers  of  scatterers  to  get  to  the  back  layer  (each  layer 
can  be  seen  independently  in  Figure  4.22).  Port  2  excitation  means  that  the  back  layer  is 
the  first  layer  experienced  by  the  signal.  Since  the  EM  held  has  not  undergone  previous 
scattering  events,  the  magnitude  of  localization  is  lower  for  this  propagation  direction. 
Using  the  same  reasoning  as  for  the  back  layer,  for  the  isolated  cubes  in  Figure  4.25 
(closest  to  Port  1  in  Figure  4.23),  Port  2  excitation  yields  higher  localization  because  the 
signal  travels  through  more  scatterers  before  reaching  the  isolated  cubes. 
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Also  from  Figure  4.25,  it  is  seen  that  the  lowest  overall  localization  magnitude  occurs 
when  the  isolated  cubes  (i.e.  few  cubes  in  a  layer)  are  the  first  scatterers  encountered  by 
the  EM  pnlse.  These  are  the  isolated  cubes  closest  to  Port  1  in  Figure  4.23  with  Port 
1  excitation.  Conversely,  the  highest  overall  localization  magnitude  occurs  for  the  back 
cubes  closest  to  Port  2  in  Figure  4.23  with  Port  1  excitation  (at  the  location  encompassed 
by  the  circle  in  Figure  4.24).  In  this  situation  there  are  the  most  scatterers  within  the 
layer  and  the  signal  has  traveled  through  the  highest  number  of  scatterers  to  reach  this 
back  point  close  to  Port  2. 

In  [56],  energy  localization  was  found  on  the  edges  and  corners  of  a  single  irregular 
inclusion  crystal.  This  chapter  studied  energy  density  behavior  for  mixtures  with  many 
disordered  irregular,  as  well  as  organized  cube,  inclusions.  As  the  EM  signal  travels 
through  more  scatterers,  the  peak  energy  density  has  been  shown  to  increase  even  if 
the  inclusions  are  small  compared  to  wavelength.  These  results  suggest  that  a  direct 
relationship  exists  between  the  disorder  introduced  by  an  increasing  number  of  inclusions 
and  localization  magnitude. 

4.8  Summary 

The  effective  EM  properties  for  mixtures  involving  many  disordered  crystals  were  stud¬ 
ied.  In  particular,  the  crystal  shape  was  found  to  have  no  appreciable  impact  on  ef¬ 
fective  permittivity  or  energy  localization  behavior.  The  effective  permittivity  derived 
from  simulations  deviates  from  that  calculated  using  conventional  Maxwell  Garnett  and 
Bruggeman  mixing  theories.  This  departure  is  especially  significant  when  there  is  a  high 
permittivity  contrast  between  the  inclusions  and  the  embedding  matrix.  This  is  so  when 
the  inclusions  have  irregular  or  regular  shapes.  For  mixtures  with  hundreds  of  crystals, 
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the  effect  of  individual  inclusion  shape  and  orientation  on  effective  permittivity  appar¬ 
ently  averages  out.  As  expected,  inclusions  in  the  form  of  a  plate  yield  either  maximum 
or  minimum  effective  permittivity  depending  on  the  cross-section  presented  to  the  EM 
signal.  In  general  a  strong  inverse  relationship  was  found  between  the  cross-sectional  area 
presented  by  the  inclusions  and  the  effective  permittivity  of  the  mixture.  While  estab¬ 
lished  mixing  laws  are  functions  of  individual  permittivity  values  and  volume  fraction, 
this  chapter  showed  that  cross-sectional  area  is  also  an  important  parameter  to  consider 
when  determining  the  effective  permittivity  of  mixtures  of  finite  thickness. 

Energy  was  seen  to  localize  within  mixtures  containing  many  inclusions.  The  specific 
location  of  EM  localization  is  on  edges  and  corners  of  adjacent  inclusions  that  are  close 
to  each  other  in  the  direction  of  electric  field  polarization.  Temporal  maximum  energy 
density  occurs  only  after  appreciable  scattering  has  occurred  but  before  the  EM  signal 
has  been  attenuated. 
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Chapter  5 


Summary  and  Future  Work 

5.1  Summary  of  Research 

This  dissertation  presented  an  overview  of  effective  permittivity  and  EM  localization  be¬ 
havior  for  various  types  of  composite  mixtures.  The  complexity  of  these  mixtures  was 
increased  throughout  the  dissertation.  In  Chapter  3,  the  examined  composites  included 
those  with  a  single  cube  inclusion,  a  single  irregular  inclusion,  and  many  ordered  cube 
inclusions  placed  on  a  regular  grid.  In  Chapter  4,  the  composites  were  much  more  com¬ 
plex,  and  included  mixtures  with  hundreds  of  regularly-  and  irregularly-shaped  inclusions. 
The  research  describing  the  characterization  of  composites  within  this  dissertation  can 
be  useful  for  various  applications.  These  include  current  research  areas  such  as  struc¬ 
ture  demolition  and  demining,  as  well  as  future  areas  such  as  the  study  of  Carbon  Fiber 
Reinforced  Composites  (CFRCs). 

In  Chapter  3,  two-component  composites  were  examined  consisting  of  cube  and 
irregularly-shaped  inclusions.  Time-domain  EM  modeling  was  used  to  calculate  effective 
permittivities  and  to  examine  energy  localization  behavior  for  these  composites.  Both 
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high  and  low  permittivity  contrast  situations  were  used  and  the  components  studied 
included  the  explosive  material  HMX. 

With  irregularly  shaped  inclusion  crystals,  effective  permittivities  were  calculated  for 
each  of  three  orthogonal  propagation  directions  (x,  y,  and  z).  The  irregular  structures 
were  not  symmetric  with  respect  to  all  polarization  directions,  and  so  the  effective  permit¬ 
tivity  showed  directional  dependency.  However,  for  the  same  electric  field  polarization, 
effective  permittivity  was  independent  of  forward  or  backward  propagation  direction.  To 
calculate  effective  permittivity,  the  composites  were  placed  within  a  TEM  environment. 
This  environment  was  created  by  setting  one  pair  of  opposite  faces  as  electric  walls,  an¬ 
other  as  magnetic  walls,  and  the  last  as  open  boundaries  that  eliminated  reflections  back 
into  the  TEM  environment. 

The  first  examined  situation  involved  HMX  inclusions,  representing  low  permittivity 
contrast  between  inclusion  and  matrix  material.  In  situations  of  low  permittivity  contrast, 
the  shape  of  the  inclusion  had  little  impact  on  effective  permittivity  up  to  a  volume 
fraction  of  0.45.  With  this  low  permittivity  contrast,  the  Maxwell  Garnett  and  Bruggeman 
mixing  theories  could  be  used  up  to  a  volume  fraction  of  0.45  with  1%  error.  With  high 
permittivity  contrast,  significant  differences  were  found  from  the  classical  mixing  rules, 
Maxwell  Garnett  and  Bruggeman,  for  volume  fractions  above  0.05.  So,  for  high  filling 
factors  the  traditional  mixing  theories  did  not  accurately  predict  the  effective  permittivity 
of  irregularly  shaped  mixtures  with  high  permittivity  contrast  between  inclusion  and 
surrounding  material. 

Mixtures  with  cube  inclusions  placed  on  an  ordered  grid  were  also  examined  in  Chap¬ 
ter  3.  This  highly  ordered  scenario  allowed  the  entire  range  of  volume  fractions  to  be 
simulated.  The  effective  permittivities  for  these  composites  were  found  to  be  close  to  the 
results  predicted  by  Maxwell  Garnett  for  all  volume  fractions. 
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Energy  localization  behavior  was  examined  with  negligible  differences  found  between 
mirrored  and  periodic  boundary  conditions.  For  both  types  of  boundary  conditions,  elec¬ 
tric  energy  localized  on  edges  and  corners  of  the  irregular  inclusion.  The  specific  edges 
and  corners  were  those  closest  to  the  adjacent  crystal  in  the  direction  of  electric  held 
polarization. 

Also  in  Chapter  3,  the  inclusion  crystals  were  rotated  and  the  orientation  of  the  inclu¬ 
sions  with  high  permittivity  contrast  impacted  the  effective  permittivity  of  the  composite. 
First,  a  cube  inclusion  at  10%  volume  fraction  was  rotated  around  the  x,  y,  and  z  axes 
independently.  With  the  electric  held  polarized  in  the  x  direction,  effective  permittivity 
was  determined  for  each  rotation  point.  Then,  a  20%  volume  fraction  cube  was  rotated 
around  each  axis  with  the  electric  held  also  polarized  in  the  x  direction.  Finally,  in  a 
more  complex  scenario,  rotation  was  performed  for  an  irregular  inclusion  at  20%  volume 
fraction.  The  irregular  inclusion  was  rotated  around  the  x,  y,  and  z  axes  independently 
and  the  electric  held  was  polarized  in  either  the  x,  y,  or  z  direction.  For  all  of  the  ro¬ 
tation  scenarios,  negligible  changes  in  effective  permittivity  resulted  when  the  inclusion 
was  rotated  about  the  same  axis  as  the  electric  held  polarization.  However,  variations  in 
effective  permittivity  resulted  from  rotations  about  the  other  axes. 

Plots  were  also  provided  showing  changes  in  cross-sectional  area  fraction  as  each  inclu¬ 
sion  was  rotated.  The  variation  in  the  cross-sectional  area  fraction  from  the  perspective 
of  a  given  axis  was  highly  correlated  to  changes  in  effective  permittivity  with  the  elec¬ 
tric  held  polarized  along  the  same  axis.  If  the  cross-sectional  area  fraction  changed  or 
remained  constant,  the  effective  permittivity  also  showed  the  same  behavior.  An  inverse 
relationship  between  cross-sectional  area  fraction  and  effective  permittivity  developed 
when  a  significant  portion  of  the  inclusion  crystal  approached  the  outer  box,  or  adja¬ 
cent  inclusions.  Differences  in  effective  permittivity  from  rotation  of  the  inclusion  crystal 
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correspond  to  a  probe  moving  around  a  specific  compound. 

In  Chapter  4,  the  complexity  of  the  examined  composites  was  increased  compared  to 
those  studied  in  Chapter  3.  EM  properties  of  mixtures  were  explored  further  by  examining 
structures  with  many  irregular  and  cube  inclusion  crystals  in  a  disordered  arrangement. 
The  crystals  had  both  low  and  high  permittivity  contrast  with  the  surrounding  material. 
With  low  permittivity  inclusions,  the  results  of  effective  permittivity  were  close  to  the 
predictions  of  Bruggeman  up  to  a  volume  fraction  of  40%.  For  the  same  volume  fraction, 
the  effective  permittivities  for  the  irregular  and  cube  inclusion  crystal  mixtures  differed 
by  less  than  0.5%.  So,  with  hundreds  of  disordered  inclusion  crystals  and  a  relatively  low 
permittivity  contrast,  the  effective  permittivity  of  mixtures  had  very  little  dependence 
on  inclusion  shape. 

Next,  the  crystals  were  given  a  higher  relative  permittivity.  In  this  situation  the 
effective  permittivity  results  were  within  the  Wiener  bounds  (as  expected)  and  were  also 
between  the  results  predicted  by  Maxwell  Garnett  and  Bruggeman  up  to  a  volume  fraction 
of  40%.  With  this  high  permittivity  contrast,  a  majority  of  the  results  between  irregular 
and  cube  inclusion  crystals  differed  by  1%  or  less  at  the  same  volume  fraction.  Here, 
even  with  a  high  permittivity  contrast  (28  to  1)  there  was  still  very  little  dependence  on 
inclusion  shape  for  effective  permittivity. 

The  relationship  between  overall  inclusion  structure  and  effective  permittivity  was 
also  examined.  Specifically,  investigations  were  performed  looking  into  how  the  arrange¬ 
ment  of  inclusions  impacted  the  individual  extracted  effective  permittivity  components 
corresponding  to  each  direction  of  electric  field  polarization.  The  Wiener  bounds  rep¬ 
resent  the  maximum  and  minimum  possible  effective  permittivity  for  a  composite  and 
result  from  an  inclusion  in  the  shape  of  a  plate.  The  maximum  and  minimum  possible  ef¬ 
fective  permittivities  correspond  to  capacitors  connected  in  parallel  or  series  respectively. 
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As  a  test,  an  inclusion  in  the  shape  of  a  plate  at  11%  volume  fraction  was  examined  in 
CST  Microwave  Studio.  As  expected,  the  extracted  effective  permittivities  corresponded 
to  the  results  predicted  by  the  Wiener  bounds. 

The  maximum  effective  permittivity  at  a  given  volume  fraction  resulted  when  the 
largest  face  of  the  inclusion  plate  was  parallel  to  the  electric  held  polarization  of  the 
TEM  held.  Conversely,  the  minimum  effective  permittivity  resulted  when  the  largest 
face  of  the  inclusion  plate  was  perpendicular  to  the  direction  of  electric  held  polarization. 

These  results  were  then  related  to  the  cross-sectional  areas  of  the  inclusion  shaped 
as  a  plate.  The  maximum  possible  effective  permittivity  resulted  when  the  inclusion 
had  the  minimum  possible  cross-sectional  area  normal  to  the  direction  of  the  electric 
held  polarization.  Conversely,  the  minimum  possible  effective  permittivity  occurred  when 
the  inclusion  had  the  maximum  possible  cross-sectional  area  normal  to  the  direction  of 
electric  held  polarization.  Combining  these  results,  there  was  an  inverse  relationship 
between  cross-sectional  area  of  the  inclusion  and  effective  permittivity  of  the  mixture  for 
an  inclusion  shaped  as  a  plate. 

This  examination  between  cross-sectional  area  fraction  and  effective  permittivity  was 
then  extended  to  include  mixtures  with  hundreds  of  inclusions.  These  mixtures  with 
hundreds  of  cube  or  irregularly-shaped  inclusions  had  close  groupings  of  inclusions  on 
the  bottom  of  the  outer  box,  forming  an  overall  structure  similar  to  a  plate.  The  effective 
permittivity  was  the  smallest  when  the  electric  held  was  polarized  along  the  axis  from 
which  the  inclusion  crystals  had  the  highest  cross-sectional  area  fraction.  Conversely,  the 
other  axes  had  lower  cross-sectional  area  fraction  and  higher  effective  permittivity.  Thus, 
a  strong  inverse  relationship  was  also  found  between  cross-sectional  area  fraction  and 
effective  permittivity  for  many  inclusions  with  this  plate- like  arrangement. 

Chapter  4  also  examined  EM  localization  for  two  distinct  arrangements  of  inclusions. 
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The  first  type  of  mixture  had  hundreds  of  irregularly-shaped  inclusion  crystals  of  various 
sizes.  At  a  volume  fraction  of  32.5%,  and  with  electric  field  polarized  along  the  x  direction, 
electric  energy  density  was  examined  for  forward  and  reverse  propagation  directions. 
The  electric  energy  density  localized  with  different  peak  magnitudes  depending  on  the 
propagation  direction  of  the  EM  signal.  At  specific  inclusion  edges  and  corners  close  to 
each  other  in  the  direction  of  electric  field  polarization,  the  magnitude  of  the  energy 
density  increased  as  the  EM  signal  traveled  through  more  inclusions. 

The  second  type  of  composite  had  dozens  of  cube  inclusion  crystals  of  the  same 
size,  where  the  inclusions  were  separated  into  parallel  planes  and  varying  numbers  of 
inclusions  were  placed  randomly  in  each  plane.  An  increasing  number  of  inclusions  were 
placed  between  the  outer  layers  to  gradually  increase  the  number  of  inclusions  in  each 
layer.  The  highest  energy  density  resulted  on  the  edges  of  adjacent  cubes  close  together 
after  the  EM  signal  traveled  through  the  highest  number  of  inclusions. 

In  addition  to  the  energy  density  results,  simulation  probes  enabled  the  electric  fields 
to  be  measured  throughout  the  simulation  at  various  points.  For  the  irregular  crystal 
mixture,  the  electric  field  localized  at  crystal  corners  and  edges  to  a  greater  amount  when 
the  EM  signal  previously  traveled  through  more  scatterers.  For  the  cube  crystal  mixture, 
electric  fields  were  measured  for  the  isolated  cubes  and  the  back  cubes  surrounded  by 
many  other  inclusions  (for  opposite  propagation  directions,  resulting  in  four  different 
electric  field  probe  results).  The  lowest  electric  field  localization  magnitude  among  these 
four  cases  occurred  for  the  isolated  cubes  when  these  cubes  were  the  first  encountered 
by  the  EM  pulse.  The  highest  localization  magnitude  occurred  for  the  back  cubes  where 
there  were  the  most  scatterers  within  that  layer  and  the  signal  had  traveled  through 
the  highest  number  of  scatterers  to  reach  that  back  point.  Overall,  higher  peak  electric 
field  localization  magnitude  was  reached  after  the  EM  pulse  traveled  through  a  greater 
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number  of  inclusions  for  both  irregularly-shaped  inclusions  in  a  disordered  arrangement 
and  more  organized  cube  crystals. 

An  important  conclusion  here  is  that  the  study  of  EM  localization  does  not  require 
geometrically  complex  structures,  such  as  those  created  with  the  irregular  inclusion  mix¬ 
tures  in  Chapter  4.  Simpler  arrangements  of  regularly  shaped  inclusions  exhibit  similar 
EM  hotspot  behavior,  but  with  a  significant  decrease  in  mixture  complexity. 

5.2  Future  Work 

This  dissertation  examined  various  types  of  composite  mixtures  containing  crystals  of 
differing  shapes  and  sizes.  One  area  that  can  be  expanded  upon  in  future  research  is  an 
increase  in  the  irregularity  of  the  crystals  that  are  examined.  For  example,  the  deviation 
in  crystal  shapes  from  one  crystal  to  the  next  could  be  increased. 

Complex  mixtures  with  many  irregular  crystals  were  created  with  volume  fractions 
up  to  42%.  Future  research  could  attempt  to  increase  the  maximum  volume  fractions 
simulated,  possibly  through  new  procedures  for  creating  crystalline  mixtures. 

The  research  in  this  dissertation  showed  similar  EM  localization  behavior  for  complex 
combinations  of  irregular  crystals  and  more  organized  cube  inclusions,  with  hotspots 
occurring  on  inclusion  edges  and  corners.  With  the  knowledge  that  studies  of  hotspots  can 
be  performed  with  relatively  simple  arrangements  of  inclusions,  moving  forward  further 
simplifications  could  be  made. 

One  idea  is  to  take  the  concept  of  using  EM  simulation  software  to  examine  combi¬ 
nations  of  crystals  within  a  surrounding  material,  as  explored  in  this  dissertation,  and 
make  an  abstraction  to  a  study  of  transmission  lines.  A  lossless  transmission  line  has  se¬ 
ries  inductance,  L,  represented  by  an  inductor,  and  shunt  capacitance,  C,  represented  by 
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a  capacitor.  LC  tanks  could  be  utilized  to  represent  inclusion  crystals.  Inclusion  corners 
could  be  represented  by  LCs  with  characteristic  impedance  representative  of  a  crystal,  and 
LCs  connected  together  with  this  impedance  could  represent  a  complete  crystal.  These 
‘crystals’  could  be  connected  to  each  other  by  LCs  with  the  characteristic  impedance 
of  vacuum,  representing  the  binder  material  around  the  crystals.  Randomness  could  be 
introduced  here  by  changing  the  length  of  the  transmission  lines. 

The  transmission  lines  study  can  be  expected  to  result  in  faster  computation  enabling 
such  things  as  optimum  EM  waveform  design.  Gaussian  input  pulses  were  used  in  this 
dissertation  in  order  to  study  a  range  of  frequencies.  The  behavior  of  EM  localization  in 
composites  can  be  further  understood  by  examining  different  types  of  input  waveforms, 
such  as  sinusoidal  or  rectangular  inputs,  and  focusing  on  specific  frequencies  to  search 
for  an  optimal  waveform  or  frequency  leading  to  high  localization.  The  study  of  energetic 
materials  in  particular  would  be  greatly  aided  with  research  showing  increased  hotspot 
magnitude  due  to  a  certain  type  of  excitation  signal. 
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Appendix  A 


MATLAB  Code 

A.l  Code  to  Calculate  Effective  Permittivity 

This  section  contains  an  example  of  the  MATLAB  code  used  to  import  the  simulated 
S'-parameters  from  CST  Microwave  Studio  and  then  calculate  effective  permittivity.  This 
code  is  used  for  the  data  in  Figure  3.9. 

clear  all; 

load  HMXTwentyEyRl.txt 
load  HMXTwentyEyTT.txt 

Sll_y_20  =  HMXTwentyEyRl ; 

S21_y_20  =  HMXTwentyEyTl; 
r_y_20  =  zeros (1001 , 1) ; 
t_y_20  =  zeros (1001 , 1) ; 
n_y_20  =  zeros (1001 , 1) ; 
z_y_20  =  zeros (1001 , 1)  ; 
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for  i  =  1:1001 


r_y_20(i)  =  Sll_y_20 (i , 2) . *exp (1 j *Sll_y_20 (i , 3) * (pi/180) ) ; 

end 

for  i  =  1:1001 

t_y_20 (i)  =  S21_y_20(i,2) . *exp (lj *S21_y_20 (i , 3) * (pi/180) ) ; 

end 

for  i  =  1:1001 

f _y_20 (i)  =  S21_y_20 (i , 1) ; 

end 

load  HMXTwentyEyR2Fixed.txt 
load  HMXTwentyEyT2Fixed.txt 

Sll_y2_20  =  HMXTwentyEyR2Fixed; 

S21_y2_20  =  HMXTwentyEyT2Fixed; 
r_y2_20  =  zeros (1001 , 1) ; 
t_y2_20  =  zeros (1001 , 1) ; 
n_y2_20  =  zeros (1001 , 1) ; 
z_y2_20  =  zeros (1001 , 1) ; 

for  i  =  1:1001 

r_y2_20(i)  =  Sll_y2_20 (i , 2) . *exp (1 j *Sll_y2_20 (i , 3) * (pi/180) ) 
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end 


for  i  =  1:1001 

t_y2_20 (i)  =  S21_y2_20 (i , 2) . *exp (1 j  *S21_y2_20 (i , 3) * (pi/ 180) ) 

end 

for  i  =  1:1001 

f _y2_20 (i)  =  S21_y2_20 (i , 1) ; 

end 

load  HMXTwentyExRlFixed.txt 
load  HMXTwentyExTlFixed.txt 

Sll_x_20  =  HMXTwentyExRlFixed; 

S21_x_20  =  HMXTwentyExTIFixed; 
r_x_20  =  zeros (1001 , 1) ; 
t_x_20  =  zeros (1001 , 1) ; 
n_x_20  =  zeros (1001 , 1) ; 
z_x_20  =  zeros (1001 , 1) ; 

for  i  =  1:1001 

r_x_20(i)  =  Sll_x_20(i ,2) . *exp(l j*Sll_x_20(i , 3) * (pi/ 180) ) ; 

end 

for  i  =  1:1001 
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t_x_20(i)  =  S21_x_20 (i , 2) . *exp (1 j  *S21_x_20 (i , 3) * (pi/ 180) ) ; 

end 

for  i  =  1:1001 

f _x_20(i)  =  S21_x_20(i , 1) ; 

end 

load  HMXTwentyExR2Fixed.txt 
load  HMXTwentyExT2Fixed.txt 

Sll_x2_20  =  HMXTwentyExR2Fixed; 

S21_x2_20  =  HMXTwentyExT2Fixed; 
r_x2_20  =  zeros (1001 , 1) ; 
t_x2_20  =  zeros (1001 , 1) ; 
n_x2_20  =  zeros (1001 , 1) ; 
z_x2_20  =  zeros (1001 , 1) ; 

for  i  =  1:1001 

r_x2_20(i)  =  Sll_x2_20 (i , 2) . *exp (lj *Sll_x2_20 (i , 3) * (pi/180) ) ; 

end 

for  i  =  1:1001 

t_x2_20(i)  =  S21_x2_20 (i , 2) . *exp (lj *S21_x2_20 (i , 3) * (pi/180) ) ; 

end 
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for  i  =  1:1001 


f _x2_20(i)  =  S21_x2_20 (i , 1) ; 

end 


load  HMXTwentyEzRlFixed.txt 
load  HMXTwentyEzTlFixed.txt 

Sll_z_20  =  HMXTwentyEzRlFixed; 

S21_z_20  =  HMXTwentyEzTIFixed; 
r_z_20  =  zeros (1001 , 1) ; 
t_z_20  =  zeros (1001 , 1) ; 
n_z_20  =  zeros (1001 , 1) ; 
z_z_20  =  zeros (1001 , 1) ; 

for  i  =  1:1001 

r_z_20(i)  =  Sll_z_20(i ,2) . *exp(l j*Sll_z_20(i , 3) * (pi/ 180) ) ; 

end 

for  i  =  1:1001 

t_z_20(i)  =  S21_z_20 (i , 2) . *exp (1 j  *S21_z_20 (i , 3) * (pi/ 180) ) ; 

end 

for  i  =  1:1001 

f _z_20(i)  =  S21_z_20(i , 1) ; 
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end 


load  HMXTwentyEzR2Fixed.txt 
load  HMXTwentyEzT2Fixed.txt 

Sll_z2_20  =  HMXTwentyEzR2Fixed; 

S21_z2_20  =  HMXTwentyEzT2Fixed; 
r_z2_20  =  zeros (1001 , 1) ; 
t_z2_20  =  zeros (1001 , 1) ; 
n_z2_20  =  zeros (1001 , 1) ; 
z_z2_20  =  zeros (1001 , 1) ; 

for  i  =  1:1001 

r_z2_20(i)  =  Sll_z2_20 (i , 2) . *exp (lj *Sll_z2_20 (i , 3) * (pi/180) ) ; 

end 

for  i  =  1:1001 

t_z2_20(i)  =  S21_z2_20 (i , 2) . *exp (lj *S21_z2_20 (i , 3) * (pi/180) ) ; 

end 

for  i  =  1:1001 

f _z2_20(i)  =  S21_z2_20  (i ,  1) ; 

end 

c=299 .8; 
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d_Ex_20  =  0.260; 


d_Ey_20  =  0.260; 
d_Ez_20  =  0.260; 

for  i  =  1:1001 

n_y_20(i)  =  (acos( (l/(2*t_y_20(i) ) )*(l-( (r_y_20(i) ) . “2  - 
(t_y_20(i) ) . ~2) ) )+(2*pi* (0) ) ) . / ( (2*pi*f _y_20(i) *d_Ey_20) . /c) ; 
z_y_20 (i)  =  sqrt ( ( (1+ (r_y_20 (i) ) ) .~2  - 

(t_y_20 (i) ) . ~2) . / ( (1- (r_y_20(i) ) ) . ~2  -  (t_y_20(i) ) . ~2) ) ; 

end 


for  i  =  1:1001 

n_y2_20(i)  =  (acos ( (1/ (2*t_y2_20 (i) ) ) * (1- ( (r_y2_20 (i) ) . ~2  - 
(t_y2_20(i)) .“2)))+(2*pi*(0))) ./((2*pi*f_y2_20(i)*d_Ey_20) ,/c) 
z_y2_20 (i)  =  sqrt ( ( (1+ (r_y2_20 (i) )  )  .  ~2  - 

(t_y2_20 (i) ).~2)./((l- (r_y2_20 (i) )  )  .  ~2  -  (t_y2_20 (i) ) . ~2) ) ; 

end 


for  i  =  1:1001 

n_x_20(i)  =  (acos((l/(2*t_x_20(i)))*(l-((r_x_20(i)).“2  - 
(t_x_20 (i) ) . ~2) ) )+(2*pi* (0) ) ) . / ( (2*pi*f _x_20 (i) *d_Ex_20) . /c) ; 
z_x_20(i)  =  sqrt ( ( (1+ (r_x_20 (i) ) ) .~2  - 

(t_x_20(i)) . ~2) ,/((l-(r_x_20(i))) . ~2  -  (t_x_20(i) ) . "2) ) ; 

end 
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for  i  =  1:1001 


n_x2_20(i)  =  (acos((l/ (2*t_x2_20(i)))*(l-((r_x2_20(i)) . ~2  - 
(t_x2_20(i)) .~2)))+(2*pi*(0))) . / ((2*pi*f_x2_20(i)*d_Ex_20) ,/c) ; 
z_x2_20(i)  =  sqrt(((l+(r_x2_20(i))) .“2  - 

(t_x2_20(i) ) .  ~2) . / ( (l-(r_x2_20(i) ) ) . ~2  -  (t_x2_20(i) ) . ~2) ) ; 

end 


for  i  =  1:1001 

n_z_20(i)  =  (acos((l/(2*t_z_20(i)))*(l-((r_z_20(i)) . "2  - 
(t_z_20(i) ) . ~2) ) )+(2*pi* (0) ) ) . / ( (2*pi*f _z_20(i) *d_Ez_20) ,/c) ; 
z_z_20(i)  =  sqrt (( (1+ (r_z_20 (i) ) ) . ~2  - 

(t_z_20(i) ) . ~2) . /( (1- (r_z_20(i) ) ) . ~2  -  (t_z_20 (i) ) . ~2) ) ; 

end 


for  i  =  1:1001 

n_z2_20(i)  =  (acos((l/(2*t_z2_20(i)))*(l-((r_z2_20(i)) . “2  - 
(t_z2_20(i)) .“2)))+(2*pi*(0))) ./ ( (2*pi*f _z2_20 (i) *d_Ez_20) ,/c) ; 
z_z2_20(i)  =  sqrt(((l+(r_z2_20(i))) . ~2  - 

(t_z2_20(i) ) . ~2) . / ( (l-(r_z2_20(i) ) ) . ~2  -  (t_z2_20(i) ) . ~2) ) ; 

end 

•/.Plot  effective  permittivity  using  only  positive  propagation  directions 
plot (f _z_20 , ( 1/3) * (real (n_y_20 . /z_y_20) +real (n_x2_20 . /z_x2_20) 

+real (n_z_20 . /z_z_20) ) , Jb’) ; 
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axis ([0.1  50  0  3] ) ; 
xlabel ( ’ f  (GHz) ’ ) ; 

ylabel(’Effective  relative  permittivity ’ ) ; 

title(’Real  effective  relative  part  of  epsilon,  Filling  factor  =  20%, 
No  rotation,  Positive  Excited  Directions ’ ) ; 
grid  on 

•/.Plot  real  part  of  effective  permittivity  using  all  6  (positive  and 
negative)  propagation  directions 
figure (2) 

plot (f _z_20 , (1/6) * (real (n_y_20 . /z_y_20) +real (n_y2_20 . /z_y2_20) + 
real (n_x_20 . /z_x_20)+real (n_x2_20 . /z_x2_20)+real (n_z_20 . /z_z_20) 

+real (n_z2_20 . /z_z2_20) ) , ’b’ ) ; 
axis ([0.1  50  0  3] ) ; 
grid  on 

xlabel ( ’ f  (GHz) ’ ) ; 

ylabel(’Real  effective  Relative  Permittivity ’ ) ; 

title(’Real  effective  relative  part  of  epsilon,  Filling  Factor  =  20%, 
No  rotation,  6  Excited  Directions’); 
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A. 2  Code  to  Plot  Effective  Permittivity  and  Mixing 
Theories 

This  example  MATLAB  code  is  used  to  plot  the  results  of  effective  permittivity  versus 
the  Maxwell  Garnett  and  Bruggeman  mixing  theories  as  seen  in  Figure  3.10. 

clear  all; 
f  =0:0.001:0.6; 
el  =  28-1 j  *0 . 0252 ; 
e2  =  l-lj*0; 

e_eff  =(e2  +(3 . *f . *e2 . *(el-e2) ) . /(el+2 . *e2-f . * (el-e2) )) ; 

plot(f ,real(e_eff ) , ’-k’ ) ; 
xlabel ( ’ Volume  fraction’) 

ylabel(’Real  effective  relative  permittivity’) 
grid  on 
hold  on 
a  =  2; 

b  =  (el-3 . *el . *f-2 . *e2+3 . *e2 . *f ) 
c  =  (-e2 . *el) ; 

e_ef f 2  =  (-b  +  sqrt(b.~2  -  4.*a.*c)) ./(2.*a) ; 
plot(f ,real(e_eff2) , ’ — k’ ) 
hold  on 

x_zirconia_pos  =  [0.1  0.2  0.25  0.3  0.35  0.4  0.45]; 
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y_zirconia_pos  =  [1.355  1.801  2.06  2.377  2.731  3.193  3.807]; 
x_zirconia_all_six  =  [0.1  0.2  0.25  0.3  0.35  0.4  0.45]; 
y_zirconia_all_six  =  [1.356  1.801  2.059  2.379  2.73  3.193  3.805]; 

plot (x_zirconia_pos ,y_zirconia_pos , ’ sk’ ) ; 
axis (  [0  0.47  0  12]); 
hold  on 

plot (x_zirconia_all_six,y_zirconia_all_six, ’+k’ ) ; 

legend ( ’ Maxwell  Garnett ’, ’Bruggeman’ , ’Finite  Integration  Technique  (FIT), 
positive  propagation  directions’ , ’Finite  Integration  Technique  (FIT), 
six  propagation  directions’); 
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